Challenging assumptions: meta-analysis reveals dispersal costs are not universal.

Supplmentary Material

Author

Martinig, A. R., S. L. P. Burk, Y. Yang, M. Lagisz, and S. Nakagawa

Published

April 23, 2025

Things to do for April

  1. Rerun all the models
  2. Add some text so the reader can understand
  3. Go through code and add more annotations
  4. Ask questions about things you do not understand

1 Setting up

1.1 Load packages

Code
#force install pacman and orchaRd

#install.packages("pacman")
#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)
#devtools::install_github("daniel1noble/orchaRd", force = TRUE)

#load all other packages
pacman::p_load(tidyverse, # tidy family and related pacakges below
               kableExtra, # nice tables
               MuMIn,  # multi-model inference
               emmeans, # post-hoc tests
               gridExtra, # may not use this
               pander,   # nice tables
               metafor,  # package for meta-analysis
               ape,      # phylogenetic analysis
               ggtree,
               ggstance,
               MuMIn,  # multi-model inference
               patchwork,   # putting ggplots together - you need to install via devtool
               here,         # making reading files easy
               orchaRd, # plotting orchard plots
               matrixcalc, # matrix calculations
               cowplot #multipanel plots
)

eval(metafor:::.MuMIn)

1.2 Load data: paper and tree data

Code
#the dataset with effect sizes
dat <- read.csv(here("data", "clean_data.csv"))

#dim(dat)
#head(dat)

#the phylogenetic tree
tree <- read.tree(here("data", "species_tree.tre"))

1.3 Calculate effect sizes and sampling variances using custom functions

Code
# function to calculate effect sizes
# Zr - correlation
# there is always n

effect_size <- function(m1, m2, sd1, sd2, n1, n2, n, # 12 arguments
                        est , se, p_val, direction, method){

  if(method == "mean_method"){
  
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s_pool <- sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r    
    
  }else if(method == "count_method"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- sqrt(m1)
    s2 <- sqrt(m2)
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r     
    
  }else if(method == "percent_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "percent_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    sd1 <- sd1/100
    sd2 <- sd2/100
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1))
    m2 <- asin(sqrt(m2))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "t_method1"){
    

    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "t_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "F_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
  
  }else if(method == "F_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "correlation_method1"){
    
    r <- est
    
  }else if(method == "correlation_method2"){
    
    r <- 2*sin((pi/6)*est)
    
  }else if(method == "estimate_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  }else if(method == "estimate_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  } 
  
    if(r >= 1){
    # if over 1, we use 0.99
    Zr <- atanh(0.99)

    }else if(r <= -1){

    Zr <- atanh(-0.99) # r = correlation

    } else {

    Zr <- atanh(r) # r = correlation

    }
  
  VZr <- 1 /(n - 3)
  
  data.frame(ri = r, yi = Zr , vi = VZr)
  
}


##########

#' Title: Contrast name generator
#'
#' @param name: a vector of character strings
cont_gen <- function(name) {
  combination <- combn(name, 2)
  name_dat <- t(combination)
  names <- paste(name_dat[, 1], name_dat[, 2], sep = "-")
  return(names)
}

#' @title get_pred1: intercept-less model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred1 <- function (model, mod = " ") {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, "")))
  len <- length(name)
  
   if (len != 1) {
        newdata <- diag(len)
        pred <- metafor::predict.rma(model, 
                                     newmods = newdata,
                                     tau2.levels = 1:len)
    }
    else {
        pred <- metafor::predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title get_pred2: normal model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred2 <- function (model, mod = " ") {
  name <- as.factor(str_replace(row.names(model$beta), 
                                paste0("relevel", "\\(", mod,", ref = name","\\)"),""))
  len <- length(name)
  
  if(len != 1){
  newdata <- diag(len)
  pred <- predict.rma(model, intercept = FALSE, newmods = newdata[ ,-1])
  }
  else {
    pred <- predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title mr_results
#' @description Function to put results of meta-regression and its contrasts
#' @param res1: data frame 1
#' @param res1: data frame 2
mr_results <- function(res1, res2) {
  restuls <-tibble(
    `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)),
    Estimate = c(res1$estimate, res2$estimate),
    `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL),
    `Upper CI  [0.975]` = c(res1$upperCL, res2$upperCL),
    `P value` = c(res1$pval, res2$pval),
    `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR),
    `Upper PI  [0.975]` = c(res1$upperPR, res2$upperPR),
  )
}


#' @title all_models
#' @description Function to take all possible models and get their results
#' @param model: intercept-less model
#' @param mod: the name of a moderator 

all_models <- function(model, mod = " ", type = "normal") {
  
  # getting the level names out
  level_names <- levels(factor(model$data[[mod]]))
  dat2 <- model$data
  mod <- mod

  # creating a function to run models
  run_rma1 <- function(name) {
      
    # variance covarinace matrix for sampling errors
    VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)
    #VCV1 <- nearPD(VCV1)$mat
      
    rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }

    run_rma2 <- function(name) {
    
            # variance covarinace matrix for sampling errors
            VCVa <- vcalc(vi = dat2$vi_abs, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)
            #VCVa <- nearPD(VCVa)$mat
               
               rma.mv(abs_yi2, V = VCVa,
               mods = ~relevel(dat2[[mod]], ref = name),
                      random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
    }
    
    run_rma3 <- function(name) {
      
              # variance covarinace matrix for sampling errors
                VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, obs = dat$effectID, rho = 0.5)
                #VCV1 <- nearPD(VCV1)$mat
               
                rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list( # hetero in relation to mod
                            formula(paste("~", mod, "| effectID")),
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     struct = "DIAG",
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }


# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])
# use different specifications for aboslute and hetero models 
if (type == "normal"){

    model_all <- purrr::map(level_names[-length(level_names)], run_rma1)

  }else if(type == "absolute"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma2)
    
  }else if(type == "hetero"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma3)
    
  }
  
  # getting estimates from intercept-less models (means for all the groups)
  res1 <- get_pred1(model, mod = mod)
  
  # getting estiamtes from all contrast models
  res2_pre <- purrr::map(model_all, ~ get_pred2(.x, mod = mod))
  
  # a list of the numbers to take out unnecessary contrasts
  contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1))
  res2 <- purrr::map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) 
  # creating a table
  res_tab <- mr_results(res1, res2) %>% 
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")
  
  # results
  res_tab

}

# absolute value analyses

# folded mean
folded_mu <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_mu
} 

# folded variance
folded_v <-function(mean, variance){
  mu <- mean
  sigma <- sqrt(variance)
  fold_mu <- sigma*sqrt(2/pi)*exp((-mu^2)/(2*sigma^2)) + mu*(1 - 2*pnorm(-mu/sigma))
  fold_se <- sqrt(mu^2 + sigma^2 - fold_mu^2)
  # adding se to make bigger mean
  fold_v <-fold_se^2
  fold_v
} 

1.4 Preparing final dataset

2 Intercept meta-analytic model

2.1 Phylogenetic multilevel models

Code
# takes a while to run
mod <- rma.mv(yi = yi, 
              V = VCV,
              mod = ~ 1,
              data = dat,
              random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned,
                            ~ 1 | phylogeny),
              R= list(phylogeny = cor_tree),
              test = "t",
              sparse = TRUE)

# saving the runs
saveRDS(mod, here("Rdata", "mod.RDS"))
Code
# getting the saved model
mod <- readRDS(here("Rdata", "mod.RDS"))

summary(mod)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.0913   558.1826   568.1826   590.7488   568.2724   

Variance Components:

            estim    sqrt  nlvls  fixed           factor    R 
sigma^2.1  0.1105  0.3324    675     no         effectID   no 
sigma^2.2  0.0110  0.1050    202     no          paperID   no 
sigma^2.3  0.0287  0.1694    146     no  species_cleaned   no 
sigma^2.4  0.0000  0.0000    146     no        phylogeny  yes 

Test for Heterogeneity:
Q(df = 674) = 390155258.0414, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0301  0.0303  -0.9947  674  0.3202  -0.0896  0.0294    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# Phylogeney accounts for nothing (0.00%), so we can take it out. 
round(i2_ml(mod),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.74              71.90               7.17              18.67 
      I2_phylogeny 
              0.00 
Code
pred_mod <- predict.rma(mod)

  estimate <- pred_mod$pred
  lowerCL <- pred_mod$ci.lb
  upperCL <- pred_mod$ci.ub 
  lowerPR <- pred_mod$cr.lb
  upperPR <- pred_mod$cr.ub 
  
  table_mod <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod <- table_mod %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_mod as RDS
saveRDS(tab_mod, here("Rdata", "tab_mod.RDS"))

Table S1. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for intercept meta-analytic model with all five random effects.

Code
tab_mod <-readRDS(here("Rdata", "tab_mod.RDS"))

tab_mod
Fixed effect estimate lowerCL upperCL pval lowerPR upperPR
intercept -0.03 -0.09 0.029 0.32 -0.793 0.733

We remove phylogeny as a random effect from our meta-analytic model because phylogeny accounts for nothing (0%)

Code
mod1 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)

summary(mod1)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.2890   486.5780   494.5780   512.7536   494.6359   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2732    696     no         effectID 
sigma^2.2  0.0212  0.1457    206     no          paperID 
sigma^2.3  0.0179  0.1337    148     no  species_cleaned 

Test for Heterogeneity:
Q(df = 695) = 15636.5654, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0240  0.0299  -0.8029  695  0.4223  -0.0827  0.0347    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(i2_ml(mod1),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.06              63.69              18.11              15.26 
Code
reduced_intercept<-orchard_plot(mod1, xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)
Code
pred_mod1 <- predict.rma(mod1)

  estimate <- pred_mod1$pred
  lowerCL <- pred_mod1$ci.lb
  upperCL <- pred_mod1$ci.ub 
  lowerPR <- pred_mod1$cr.lb
  upperPR <- pred_mod1$cr.ub 
  
  table_mod1 <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod1$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod1 <- table_mod1 %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_mod1 as RDS
saveRDS(tab_mod1, here("Rdata", "tab_mod1.RDS"))

Table S2. Mean estimate with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval) for reduced intercept meta-analytic model without phylogeny as a random effect.

Code
tab_mod1 <-readRDS(here("Rdata", "tab_mod1.RDS"))

tab_mod1
Fixed effect estimate lowerCL upperCL pval lowerPR upperPR
intercept -0.024 -0.083 0.035 0.422 -0.689 0.641

3 Uni-moderator meta-regression models

For each of our moderators, we run a uni-moderator meta-regression model

Code
mod_class <- rma.mv(yi = yi, V = VCV,
               mod = ~ species_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_class)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-240.9746   481.9491   499.9491   540.7793   500.2138   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0737  0.2715    696     no         effectID 
sigma^2.2  0.0154  0.1241    206     no          paperID 
sigma^2.3  0.0302  0.1738    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 690) = 15096.1494, p-val < .0001

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 690) = 0.5620, p-val = 0.7607

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
species_classActinopterygii    0.0239  0.1665   0.1434  690  0.8860  -0.3030 
species_classArachnida         0.0426  0.2711   0.1573  690  0.8751  -0.4897 
species_classAves             -0.0262  0.0361  -0.7254  690  0.4684  -0.0970 
species_classInsecta           0.0971  0.1395   0.6960  690  0.4866  -0.1768 
species_classMammalia         -0.0544  0.0451  -1.2067  690  0.2280  -0.1430 
species_classReptilia          0.1153  0.1440   0.8007  690  0.4236  -0.1674 
                              ci.ub    
species_classActinopterygii  0.3508    
species_classArachnida       0.5749    
species_classAves            0.0447    
species_classInsecta         0.3709    
species_classMammalia        0.0341    
species_classReptilia        0.3980    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_class)
   R2_marginal R2_conditional 
    0.00899031     0.38789761 
Code
orchard_plot(mod_class, mod = "species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S2. Mean fitness effect of different taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_class <- all_models(mod_class,  mod = "species_class", type = "normal")

# saving tab_class as RDS
saveRDS(tab_class, here("Rdata", "tab_class.RDS"))

Table S3. Mean estimates of different taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_class <-readRDS(here("Rdata", "tab_class.RDS"))

tab_class
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Actinopterygii 0.001 -0.344 0.346 0.996 -0.850 0.851
Arachnida 0.039 -0.531 0.610 0.892 -0.925 1.004
Aves -0.027 -0.101 0.046 0.464 -0.808 0.753
Insecta 0.094 -0.202 0.391 0.533 -0.738 0.926
Mammalia -0.071 -0.163 0.022 0.133 -0.853 0.712
Reptilia 0.067 -0.219 0.353 0.645 -0.761 0.895
Actinopterygii-Arachnida 0.038 -0.626 0.703 0.910 -0.984 1.061
Actinopterygii-Aves -0.028 -0.378 0.321 0.873 -0.881 0.824
Actinopterygii-Insecta 0.093 -0.358 0.545 0.685 -0.806 0.993
Actinopterygii-Mammalia -0.071 -0.423 0.280 0.690 -0.925 0.782
Actinopterygii-Reptilia 0.066 -0.376 0.508 0.769 -0.828 0.961
Arachnida-Aves -0.067 -0.640 0.507 0.819 -1.033 0.899
Arachnida-Insecta 0.055 -0.586 0.696 0.866 -0.953 1.063
Arachnida-Mammalia -0.110 -0.685 0.465 0.708 -1.077 0.857
Arachnida-Reptilia 0.028 -0.607 0.663 0.931 -0.976 1.032
Aves-Insecta 0.122 -0.180 0.424 0.429 -0.712 0.956
Aves-Mammalia -0.043 -0.147 0.061 0.417 -0.827 0.741
Aves-Reptilia 0.095 -0.194 0.384 0.520 -0.735 0.924
Insecta-Mammalia -0.165 -0.470 0.140 0.289 -1.000 0.670
Insecta-Reptilia -0.027 -0.434 0.380 0.896 -0.905 0.850
Mammalia-Reptilia 0.138 -0.152 0.428 0.352 -0.692 0.967
Code
mod_type <- rma.mv(yi = yi,  
               V = VCV,
               mod = ~ study_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_type)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.0019   486.0038   496.0038   518.7161   496.0910   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0745  0.2729    696     no         effectID 
sigma^2.2  0.0205  0.1430    206     no          paperID 
sigma^2.3  0.0198  0.1406    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 15590.0770, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 0.3754, p-val = 0.6872

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
study_typenatural        -0.0233  0.0306  -0.7623  694  0.4461  -0.0834  0.0368 
study_typesemi-natural   -0.0412  0.0668  -0.6171  694  0.5374  -0.1723  0.0899 
                          
study_typenatural         
study_typesemi-natural    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_type)
   R2_marginal R2_conditional 
  0.0002521658   0.3509052501 
Code
orchard_plot(mod_type, mod = "study_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S3. Mean fitness effect of natural and semi-natural study types. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_type <- all_models(mod_type,  mod = "study_type", type = "normal")

# saving tab_type as RDS
saveRDS(tab_type, here("Rdata", "tab_type.RDS"))

Table S4. Mean estimates of natural and semi-natural study types with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_type <-readRDS(here("Rdata", "tab_type.RDS"))

tab_type
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Natural -0.028 -0.089 0.033 0.364 -0.795 0.738
Semi-natural -0.062 -0.196 0.072 0.366 -0.838 0.714
Natural-Semi-natural -0.034 -0.165 0.098 0.616 -0.809 0.742
Code
mod_design <- rma.mv(yi = yi, V = VCV,
               mod = ~ study_design - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_design)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.9601   485.9202   495.9202   518.6325   496.0074   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2730    696     no         effectID 
sigma^2.2  0.0233  0.1526    206     no          paperID 
sigma^2.3  0.0158  0.1256    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 15620.6163, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 0.4424, p-val = 0.6427

Model Results:

                      estimate      se     tval   df    pval    ci.lb   ci.ub 
study_designcontrast   -0.0189  0.0310  -0.6084  694  0.5431  -0.0798  0.0420 
study_designgradient   -0.0415  0.0459  -0.9043  694  0.3662  -0.1315  0.0486 
                        
study_designcontrast    
study_designgradient    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_design)
   R2_marginal R2_conditional 
  0.0008110263   0.3444046941 
Code
orchard_plot(mod_design, mod = "study_design", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S4. Mean fitness effect of gradient and contrast study designs. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_design <- all_models(mod_design,  mod = "study_design", type = "normal")

# saving tab_design as RDS
saveRDS(tab_design, here("Rdata", "tab_design.RDS"))

Table S5. Mean estimates of gradient and contrast study designs with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_design <-readRDS(here("Rdata", "tab_design.RDS"))

tab_design
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Contrast 0.007 -0.033 0.047 0.735 -0.575 0.589
Gradient -0.039 -0.106 0.028 0.257 -0.624 0.546
Contrast-Gradient -0.025 -0.107 0.057 0.553 -0.758 0.709
Code
mod_disp <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ dispersal_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_disp)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.7712   485.5425   497.5425   524.7886   497.6649   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2732    696     no         effectID 
sigma^2.2  0.0230  0.1516    206     no          paperID 
sigma^2.3  0.0164  0.1281    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 693) = 15475.9835, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 693) = 0.4132, p-val = 0.7436

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
dispersal_typeBoth       -0.0810  0.0804  -1.0068  693  0.3144  -0.2389  0.0769 
dispersal_typebreeding   -0.0144  0.0424  -0.3392  693  0.7346  -0.0976  0.0689 
dispersal_typenatal      -0.0218  0.0325  -0.6695  693  0.5034  -0.0855  0.0420 
                          
dispersal_typeBoth        
dispersal_typebreeding    
dispersal_typenatal       

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_disp)
   R2_marginal R2_conditional 
   0.001361912    0.346299964 
Code
orchard_plot(mod_disp, mod = "dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S5. Mean fitness effect of natal and breeding dispersal (or both dispersal types grouped). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_disp <- all_models(mod_disp,  mod = "dispersal_type", type = "normal")

# saving tab_disp as RDS
saveRDS(tab_disp, here("Rdata", "tab_disp.RDS"))

Table S6. Mean estimates of natal and breeding dispersal (or both dispersal types combined) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_disp <-readRDS(here("Rdata", "tab_disp.RDS"))

tab_disp
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.092 -0.241 0.057 0.226 -0.867 0.683
Breeding -0.013 -0.100 0.075 0.772 -0.779 0.753
Natal -0.027 -0.092 0.038 0.410 -0.791 0.736
B-Breeding 0.079 -0.084 0.242 0.343 -0.699 0.857
B-Natal 0.065 -0.087 0.216 0.403 -0.711 0.840
Breeding-Natal -0.014 -0.100 0.071 0.742 -0.780 0.751
Code
mod_phase <- rma.mv(yi = yi, V = VCV,
               mod = ~ dispersal_phase - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_phase)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.4666   484.9331   494.9331   517.6455   495.0203   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0743  0.2726    696     no         effectID 
sigma^2.2  0.0216  0.1469    206     no          paperID 
sigma^2.3  0.0180  0.1342    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 15632.6394, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 1.0905, p-val = 0.3366

Model Results:

                            estimate      se     tval   df    pval    ci.lb 
dispersal_phaseprospection   -0.0652  0.0447  -1.4585  694  0.1452  -0.1529 
dispersal_phasesettlement    -0.0141  0.0310  -0.4538  694  0.6501  -0.0750 
                             ci.ub    
dispersal_phaseprospection  0.0226    
dispersal_phasesettlement   0.0469    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_phase)
   R2_marginal R2_conditional 
   0.003449693    0.349617988 
Code
orchard_plot(mod_phase, mod = "dispersal_phase", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S6. Mean fitness effect of prospection and settlement dispersal phases. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_phase <- all_models(mod_phase,  mod = "dispersal_phase", type = "normal")

# saving tab_phase as RDS
saveRDS(tab_phase, here("Rdata", "tab_phase.RDS"))

Table S7. Mean estimates of prospection and settlement dispersal phases with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_phase <-readRDS(here("Rdata", "tab_phase.RDS"))

tab_phase
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Prospection -0.067 -0.158 0.025 0.152 -0.836 0.703
Settlement -0.022 -0.085 0.041 0.493 -0.788 0.745
Prospection-Settlement 0.045 -0.042 0.132 0.314 -0.724 0.813
Code
mod_sex <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.2724   484.5447   496.5447   523.7909   496.6672   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0741  0.2722    696     no         effectID 
sigma^2.2  0.0239  0.1547    206     no          paperID 
sigma^2.3  0.0156  0.1249    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 693) = 15546.8360, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 693) = 1.0504, p-val = 0.3696

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
sexBoth     -0.0702  0.0420  -1.6722  693  0.0949  -0.1526  0.0122  . 
sexFemale   -0.0082  0.0338  -0.2433  693  0.8078  -0.0745  0.0581    
sexMale     -0.0037  0.0349  -0.1058  693  0.9158  -0.0723  0.0649    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex)
   R2_marginal R2_conditional 
   0.005223149    0.351422447 
Code
orchard_plot(mod_sex, mod = "sex", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S7. Mean fitness effect of females, males, or both sexes grouped. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex <- all_models(mod_sex,  mod = "sex", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_sex, here("Rdata", "tab_sex.RDS"))

Table S8. Mean estimates of females, males, or both sexes grouped with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex <-readRDS(here("Rdata", "tab_sex.RDS"))

tab_sex
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.073 -0.162 0.015 0.102 -0.841 0.694
F -0.013 -0.083 0.056 0.711 -0.779 0.753
M -0.017 -0.088 0.055 0.650 -0.782 0.749
B-F 0.060 -0.032 0.152 0.199 -0.708 0.828
B-M 0.057 -0.037 0.151 0.237 -0.711 0.825
F-M -0.003 -0.062 0.055 0.908 -0.768 0.761

3.1 Higher life history grouping

Code
mod_age1 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age1)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-241.5610   483.1219   495.1219   522.3681   495.2444   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0750  0.2738    696     no         effectID 
sigma^2.2  0.0172  0.1312    206     no          paperID 
sigma^2.3  0.0219  0.1479    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 693) = 15581.3475, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 693) = 1.4496, p-val = 0.2272

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
age_class_cleanadult      -0.0190  0.0310  -0.6131  693  0.5400  -0.0799 
age_class_cleanjuvenile    0.0059  0.0513   0.1151  693  0.9084  -0.0948 
age_class_cleanmix        -0.1213  0.0607  -1.9999  693  0.0459  -0.2404 
                           ci.ub    
age_class_cleanadult      0.0419    
age_class_cleanjuvenile   0.1066    
age_class_cleanmix       -0.0022  * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age1)
   R2_marginal R2_conditional 
   0.006433797    0.346790908 
Code
orchard_plot(mod_age1, mod = "age_class_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S8. Mean fitness effect of life history stage (adult, juvenile, or mixed). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_age1 <- all_models(mod_age1,  mod = "age_class_clean", type = "normal")

# saving tab_age1 as RDS
saveRDS(tab_age1, here("Rdata", "tab_age1.RDS"))

Table S9. Mean estimates of life history stage (adult, juvenile, or mixed) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_age1 <-readRDS(here("Rdata", "tab_age1.RDS"))

tab_age1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Adult -0.004 -0.061 0.052 0.879 -0.735 0.726
Juvenile 0.012 -0.080 0.105 0.794 -0.722 0.747
Mix -0.101 -0.224 0.022 0.108 -0.839 0.638
Adult-Juvenile 0.017 -0.072 0.106 0.713 -0.717 0.751
Adult-Mix -0.096 -0.217 0.025 0.118 -0.835 0.642
Juvenile-Mix -0.113 -0.254 0.027 0.115 -0.855 0.629

3.2 Finer life stage grouping

Code
# age_class

mod_age2 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age2)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-240.3580   480.7160   500.7160   546.0684   501.0405   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0754  0.2746    696     no         effectID 
sigma^2.2  0.0183  0.1354    206     no          paperID 
sigma^2.3  0.0198  0.1407    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 689) = 15382.5080, p-val < .0001

Test of Moderators (coefficients 1:7):
F(df1 = 7, df2 = 689) = 0.9165, p-val = 0.4929

Model Results:

              estimate      se     tval   df    pval    ci.lb   ci.ub    
age_classA     -0.0112  0.0338  -0.3309  689  0.7408  -0.0775  0.0552    
age_classJ      0.0042  0.0515   0.0811  689  0.9354  -0.0969  0.1053    
age_classJA    -0.3023  0.1796  -1.6835  689  0.0927  -0.6549  0.0503  . 
age_classJY    -0.1398  0.0882  -1.5850  689  0.1134  -0.3130  0.0334    
age_classJYA   -0.0662  0.0854  -0.7745  689  0.4389  -0.2340  0.1016    
age_classY     -0.0432  0.0574  -0.7528  689  0.4519  -0.1559  0.0695    
age_classYA    -0.0375  0.0467  -0.8028  689  0.4224  -0.1292  0.0542    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age2)
   R2_marginal R2_conditional 
    0.01178388     0.34371431 
Code
orchard_plot(mod_age2, mod = "age_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S9. Mean fitness effect of life history stage at the level of precision reported (A: adult, Y: yearling, J: juvenile). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# table not created for this (not meaningful for some levels)

3.3 Higher level

Code
mod_fit1 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_higher_level - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit1)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.1277   484.2555   496.2555   523.5017   496.3779   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0744  0.2728    696     no         effectID 
sigma^2.2  0.0211  0.1453    206     no          paperID 
sigma^2.3  0.0180  0.1342    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 693) = 15506.0760, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 693) = 1.2241, p-val = 0.3000

Model Results:

                                      estimate      se     tval   df    pval 
fitness_higher_levelindirect fitness   -0.2389  0.1452  -1.6452  693  0.1004 
fitness_higher_levelreproduction       -0.0103  0.0324  -0.3193  693  0.7496 
fitness_higher_levelsurvival           -0.0338  0.0332  -1.0180  693  0.3090 
                                        ci.lb   ci.ub    
fitness_higher_levelindirect fitness  -0.5239  0.0462    
fitness_higher_levelreproduction      -0.0740  0.0533    
fitness_higher_levelsurvival          -0.0990  0.0314    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit1)
   R2_marginal R2_conditional 
   0.004159334    0.347294612 
Code
orchard_plot(mod_fit1, mod = "fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S10. Mean fitness effect of different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_fit1 <- all_models(mod_fit1,  mod = "fitness_higher_level", type = "normal")

# saving tab_fit1 as RDS
saveRDS(tab_fit1, here("Rdata", "tab_fit1.RDS"))

Table S10. Mean estimates of different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_fit1 <-readRDS(here("Rdata", "tab_fit1.RDS"))

tab_fit1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Indirect fitness -0.237 -0.572 0.098 0.166 -1.068 0.595
Reproduction -0.019 -0.084 0.047 0.580 -0.782 0.745
Survival -0.038 -0.106 0.030 0.273 -0.802 0.726
Indirect fitness-Reproduction 0.218 -0.117 0.554 0.202 -0.613 1.050
Indirect fitness-Survival 0.199 -0.136 0.534 0.244 -0.632 1.030
Reproduction-Survival -0.019 -0.079 0.040 0.525 -0.783 0.744

3.4 Finer level

Code
mod_fit2 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_metric_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit2)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-241.0975   482.1951   504.1951   554.0667   504.5856   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0740  0.2721    696     no         effectID 
sigma^2.2  0.0166  0.1287    206     no          paperID 
sigma^2.3  0.0257  0.1602    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 688) = 14978.0697, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 688) = 0.9608, p-val = 0.4656

Model Results:

                                                           estimate      se 
fitness_metric_cleanage at first reproduction               -0.0481  0.0922 
fitness_metric_cleanindirect fitness                        -0.2312  0.1459 
fitness_metric_cleanlifetime breeding success                0.0043  0.0372 
fitness_metric_cleanlifetime reproductive success           -0.0264  0.0400 
fitness_metric_cleanoffspring reproduction                   0.0341  0.0877 
fitness_metric_cleanoffspring survival                       0.0105  0.0445 
fitness_metric_cleanreproductive lifespan and/or attempts   -0.0438  0.0899 
fitness_metric_cleansurvival                                -0.0586  0.0370 
                                                              tval   df    pval 
fitness_metric_cleanage at first reproduction              -0.5218  688  0.6020 
fitness_metric_cleanindirect fitness                       -1.5847  688  0.1135 
fitness_metric_cleanlifetime breeding success               0.1166  688  0.9072 
fitness_metric_cleanlifetime reproductive success          -0.6598  688  0.5096 
fitness_metric_cleanoffspring reproduction                  0.3891  688  0.6974 
fitness_metric_cleanoffspring survival                      0.2359  688  0.8136 
fitness_metric_cleanreproductive lifespan and/or attempts  -0.4874  688  0.6262 
fitness_metric_cleansurvival                               -1.5850  688  0.1134 
                                                             ci.lb   ci.ub    
fitness_metric_cleanage at first reproduction              -0.2292  0.1329    
fitness_metric_cleanindirect fitness                       -0.5177  0.0553    
fitness_metric_cleanlifetime breeding success              -0.0688  0.0775    
fitness_metric_cleanlifetime reproductive success          -0.1050  0.0522    
fitness_metric_cleanoffspring reproduction                 -0.1381  0.2064    
fitness_metric_cleanoffspring survival                     -0.0768  0.0978    
fitness_metric_cleanreproductive lifespan and/or attempts  -0.2202  0.1326    
fitness_metric_cleansurvival                               -0.1312  0.0140    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit2)
   R2_marginal R2_conditional 
    0.01004347     0.36962120 
Code
orchard_plot(mod_fit2, mod = "fitness_metric_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S11. Mean fitness effect of different fitness components at a finer resolution. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# table not created for this (not meaningful for some levels)
Code
mod_gen <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ whose_fitness - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_gen)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.6778   485.3556   495.3556   518.0679   495.4428   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0740  0.2721    696     no         effectID 
sigma^2.2  0.0179  0.1337    206     no          paperID 
sigma^2.3  0.0240  0.1548    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 15635.4098, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 1.1603, p-val = 0.3140

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
whose_fitnessdescendant    0.0118  0.0430   0.2753  694  0.7832  -0.0725 
whose_fitnessfocal        -0.0325  0.0310  -1.0486  694  0.2947  -0.0932 
                          ci.ub    
whose_fitnessdescendant  0.0962    
whose_fitnessfocal       0.0283    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_gen)
   R2_marginal R2_conditional 
   0.002332594    0.362620116 
Code
orchard_plot(mod_gen, mod = "whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Figure S12. Mean fitness effect of focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_gen <- all_models(mod_gen,  mod = "whose_fitness", type = "normal")

# saving tab_gen as RDS
saveRDS(tab_gen, here("Rdata", "tab_gen.RDS"))

Table S11. Mean estimates of non-dispersers and dispersers and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_gen <-readRDS(here("Rdata", "tab_gen.RDS"))

tab_gen
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Individual -0.038 -0.100 0.023 0.219 -0.808 0.731
Offspring 0.005 -0.086 0.096 0.907 -0.767 0.777
Individual-Offspring 0.044 -0.035 0.123 0.275 -0.727 0.815

3.5 Normal analysis

Code
# getting absolute values for effect sizes as we expect shorter years will have more fluctuations
dat$ln_study_duration <- log(dat$study_duration+1)

mod_dur <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ ln_study_duration,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_dur)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.0283   486.0566   496.0566   518.7690   496.1438   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2731    696     no         effectID 
sigma^2.2  0.0215  0.1467    206     no          paperID 
sigma^2.3  0.0182  0.1348    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 15598.3840, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 0.0060, p-val = 0.9380

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0202  0.0599  -0.3368  694  0.7364  -0.1378  0.0974    
ln_study_duration   -0.0016  0.0211  -0.0778  694  0.9380  -0.0430  0.0398    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_dur)
   R2_marginal R2_conditional 
  2.035908e-05   3.473406e-01 
Code
bubble_plot(mod_dur,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S13. Mean fitness effect of study duration (log years). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# absolute value analysis
mod_dur_abs <- rma.mv(yi = yi_abs,
                V = VCVa,
               mod = ~ ln_study_duration,
               data = dat,
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_dur_abs)

Multivariate Meta-Analysis Model (k = 696; method: REML)

  logLik  Deviance       AIC       BIC      AICc   
 28.7611  -57.5223  -47.5223  -24.8099  -47.4351   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0326  0.1806    696     no         effectID 
sigma^2.2  0.0081  0.0900    206     no          paperID 
sigma^2.3  0.0090  0.0948    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 9232.8290, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 3.7270, p-val = 0.0539

Model Results:

                   estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt             -0.0359  0.0396  -0.9060  694  0.3653  -0.1137  0.0419    
ln_study_duration    0.0268  0.0139   1.9305  694  0.0539  -0.0005  0.0540  . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_dur_abs)
   R2_marginal R2_conditional 
    0.01232166     0.35185508 
Code
bubble_plot(mod_dur_abs,
                mod = "ln_study_duration",
                group = "paperID",
                xlab = "Study duration",
                ylab = "Effect Size: Zr",
                       g = TRUE)

Figure S14. Mean fitness effect of study duration (log years) from absolute value analysis. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).

3.6 Normal analysis

Code
mod_focus <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_focus)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.0004   486.0008   496.0008   518.7132   496.0881   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0744  0.2728    696     no         effectID 
sigma^2.2  0.0213  0.1459    206     no          paperID 
sigma^2.3  0.0188  0.1371    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 15614.7766, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 0.4197, p-val = 0.6574

Model Results:

                                 estimate      se     tval   df    pval 
fitness_main_focusNon-selective   -0.0097  0.0462  -0.2102  694  0.8336 
fitness_main_focusSelective       -0.0281  0.0313  -0.8982  694  0.3694 
                                   ci.lb   ci.ub    
fitness_main_focusNon-selective  -0.1003  0.0809    
fitness_main_focusSelective      -0.0896  0.0334    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focus)
   R2_marginal R2_conditional 
  0.0005343088   0.3503769299 
Code
orchard_plot(mod_focus, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S15. Mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_focus <- all_models(mod_focus,  mod = "fitness_main_focus", type = "normal")

# saving tab_focus as RDS
saveRDS(tab_focus, here("Rdata", "tab_focus.RDS"))

Table S14. Mean estimates of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_focus <-readRDS(here("Rdata", "tab_focus.RDS"))

tab_focus
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
N -0.019 -0.110 0.073 0.690 -0.787 0.750
Y -0.034 -0.096 0.029 0.292 -0.799 0.732
N-Y -0.015 -0.103 0.073 0.734 -0.783 0.752

3.7 Comparing normal model to heteroscedastic model

Code
# homo
mod_focus1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_focus1)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-243.8844  2860.9222   497.7688   520.4956   497.8558   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2731    696     no         effectID 
sigma^2.2  0.0214  0.1464    206     no          paperID 
sigma^2.3  0.0173  0.1314    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 15614.7766, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 694) = 0.4016, p-val = 0.6694

Model Results:

                                 estimate      se     tval   df    pval 
fitness_main_focusNon-selective   -0.0095  0.0459  -0.2066  694  0.8364 
fitness_main_focusSelective       -0.0273  0.0310  -0.8792  694  0.3796 
                                   ci.lb   ci.ub    
fitness_main_focusNon-selective  -0.0995  0.0806    
fitness_main_focusSelective      -0.0882  0.0336    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focus1)
   R2_marginal R2_conditional 
  0.0005038894   0.3419616519 
Code
# hetero
mod_focus2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ fitness_main_focus,
                 data = dat, 
                 random = list(
                   ~ fitness_main_focus | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_focus2)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-238.8875  2850.9283   489.7749   517.0470   489.8968   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0241  0.1553    206     no          paperID 
sigma^2.2  0.0163  0.1277    148     no  species_cleaned 

outer factor: effectID           (nlvls = 696)
inner factor: fitness_main_focus (nlvls = 2)

            estim    sqrt  k.lvl  fixed          level 
tau^2.1    0.0430  0.2073    164     no  Non-selective 
tau^2.2    0.0817  0.2859    532     no      Selective 

Test for Residual Heterogeneity:
QE(df = 694) = 15614.7766, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 0.2219, p-val = 0.6377

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
intrcpt                       -0.0057  0.0433  -0.1323  694  0.8948  -0.0908 
fitness_main_focusSelective   -0.0195  0.0415  -0.4711  694  0.6377  -0.1009 
                              ci.ub    
intrcpt                      0.0793    
fitness_main_focusSelective  0.0619    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
#r2_ml(mod_focusA1)

# they are not significantly different
anova(mod_focus1, mod_focus2)

        df      AIC      BIC     AICc    logLik    LRT   pval         QE 
Full     6 489.7749 517.0470 489.8968 -238.8875               15614.7766 
Reduced  5 497.7688 520.4956 497.8558 -243.8844 9.9939 0.0016 15614.7766 
Code
# homo
orchard_plot(mod_focus1, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S16. Maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# hetero
orchard_plot(mod_focus2, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S17. Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any selection bias (i.e., studies explicitly interested in the fitness effects of dispersal; Y: yes; N: no). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).

3.8 Normal analysis

Code
mod_confirm <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_confirm)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.2655   484.5309   498.5309   530.3080   498.6947   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0743  0.2726    696     no         effectID 
sigma^2.2  0.0275  0.1659    206     no          paperID 
sigma^2.3  0.0122  0.1104    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 692) = 15529.2933, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 692) = 0.3695, p-val = 0.8304

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0300  0.0829   0.3613  692  0.7180  -0.1328 
confirmation_biasCost      -0.0367  0.0389  -0.9428  692  0.3461  -0.1132 
confirmation_biasMixed     -0.0254  0.0394  -0.6460  692  0.5185  -0.1028 
confirmation_biasNeutral   -0.0008  0.0469  -0.0163  692  0.9870  -0.0929 
                           ci.ub    
confirmation_biasBenefit  0.1928    
confirmation_biasCost     0.0397    
confirmation_biasMixed    0.0519    
confirmation_biasNeutral  0.0914    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_confirm)
   R2_marginal R2_conditional 
   0.002617906    0.349759459 
Code
orchard_plot(mod_confirm, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S18. Mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_confirm <- all_models(mod_confirm,  mod = "confirmation_bias", type = "normal")

# saving tab_confirm as RDS
saveRDS(tab_confirm, here("Rdata", "tab_confirm.RDS"))

Table S15. Mean estimates of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_confirm <-readRDS(here("Rdata", "tab_confirm.RDS"))

tab_confirm
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Benefit 0.012 -0.154 0.178 0.884 -0.764 0.788
Cost -0.042 -0.122 0.037 0.296 -0.804 0.720
Mixed -0.028 -0.107 0.050 0.483 -0.790 0.734
Neutral -0.009 -0.102 0.084 0.852 -0.772 0.755
Benefit-Cost -0.055 -0.227 0.118 0.536 -0.832 0.723
Benefit-Mixed -0.040 -0.212 0.132 0.645 -0.818 0.737
Benefit-Neutral -0.021 -0.198 0.156 0.815 -0.800 0.757
Cost-Mixed 0.014 -0.081 0.110 0.770 -0.750 0.778
Cost-Neutral 0.033 -0.071 0.138 0.531 -0.732 0.799
Mixed-Neutral 0.019 -0.085 0.124 0.719 -0.746 0.784

3.9 Comparing normal model to heteroscedastic model

Code
# homo
mod_confirm1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_confirm1)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-243.5314  2860.2161   501.0628   532.8802   501.2255   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0746  0.2731    696     no         effectID 
sigma^2.2  0.0263  0.1623    206     no          paperID 
sigma^2.3  0.0110  0.1050    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 692) = 15529.2933, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 692) = 0.3600, p-val = 0.8371

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0297  0.0819   0.3623  692  0.7173  -0.1311 
confirmation_biasCost      -0.0358  0.0384  -0.9320  692  0.3516  -0.1112 
confirmation_biasMixed     -0.0241  0.0387  -0.6221  692  0.5341  -0.1001 
confirmation_biasNeutral   -0.0002  0.0463  -0.0046  692  0.9963  -0.0911 
                           ci.ub    
confirmation_biasBenefit  0.1905    
confirmation_biasCost     0.0396    
confirmation_biasMixed    0.0519    
confirmation_biasNeutral  0.0907    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_confirm1)
   R2_marginal R2_conditional 
   0.002571009    0.335545088 
Code
# hetero
mod_confirm2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ confirmation_bias,
                 data = dat, 
                 random = list(
                   ~ confirmation_bias | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_confirm2)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-224.2388  2821.6311   468.4777   513.9312   468.7989   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0288  0.1698    206     no          paperID 
sigma^2.2  0.0052  0.0719    148     no  species_cleaned 

outer factor: effectID          (nlvls = 696)
inner factor: confirmation_bias (nlvls = 4)

            estim    sqrt  k.lvl  fixed    level 
tau^2.1    0.0129  0.1135     32     no  Benefit 
tau^2.2    0.0767  0.2769    210     no     Cost 
tau^2.3    0.0565  0.2377    289     no    Mixed 
tau^2.4    0.1380  0.3715    165     no  Neutral 

Test for Residual Heterogeneity:
QE(df = 692) = 15529.2933, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 692) = 0.6362, p-val = 0.5919

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
intrcpt                     0.0503  0.0641   0.7840  692  0.4333  -0.0756 
confirmation_biasCost      -0.0921  0.0680  -1.3547  692  0.1760  -0.2255 
confirmation_biasMixed     -0.0779  0.0672  -1.1593  692  0.2467  -0.2099 
confirmation_biasNeutral   -0.0630  0.0736  -0.8552  692  0.3928  -0.2075 
                           ci.ub    
intrcpt                   0.1762    
confirmation_biasCost     0.0414    
confirmation_biasMixed    0.0540    
confirmation_biasNeutral  0.0816    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# they are not significantly different
anova(mod_confirm1, mod_confirm2)

        df      AIC      BIC     AICc    logLik     LRT   pval         QE 
Full    10 468.4777 513.9312 468.7989 -224.2388                15529.2933 
Reduced  7 501.0628 532.8802 501.2255 -243.5314 38.5851 <.0001 15529.2933 
Code
# homo
orchard_plot(mod_confirm1, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S19. Maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# hetero
orchard_plot(mod_confirm2, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Figure S20. Heteroscedastic maximum likelihood model of mean fitness effect of whether studies showed any confirmation bias (i.e., studies that stated an assumption about how dispersal affects fitness; Benefit: assumed a fitness benefit; Cost: assumed a fitness cost; Mixed: assumed fitness costs and benefits; Neutral: no assumption made). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).

4 Interaction meta-regression models (2 moderators)

Code
# the interaction between 2 categorical variables are the same as creating a new variable with 2 categories combined so we use that approach 

dat$sex_species_class <- as.factor(paste(dat$sex, dat$species_class ,sep = "_"))

mod_sex_species_class <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_species_class - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_species_class)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-236.6135   473.2270   509.2270   590.6511   510.2602   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0741  0.2722    696     no         effectID 
sigma^2.2  0.0207  0.1439    206     no          paperID 
sigma^2.3  0.0195  0.1397    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 681) = 14400.5352, p-val < .0001

Test of Moderators (coefficients 1:15):
F(df1 = 15, df2 = 681) = 0.9762, p-val = 0.4787

Model Results:

                                        estimate      se     tval   df    pval 
sex_species_classBoth_Aves               -0.0077  0.0517  -0.1480  681  0.8824 
sex_species_classBoth_Insecta            -0.1445  0.2402  -0.6015  681  0.5477 
sex_species_classBoth_Mammalia           -0.1764  0.0646  -2.7291  681  0.0065 
sex_species_classBoth_Reptilia            0.1653  0.2740   0.6034  681  0.5465 
sex_species_classFemale_Actinopterygii   -0.0592  0.1944  -0.3047  681  0.7607 
sex_species_classFemale_Arachnida         0.0410  0.2612   0.1570  681  0.8753 
sex_species_classFemale_Aves             -0.0326  0.0413  -0.7909  681  0.4293 
sex_species_classFemale_Insecta           0.2482  0.1764   1.4071  681  0.1599 
sex_species_classFemale_Mammalia         -0.0130  0.0495  -0.2632  681  0.7925 
sex_species_classFemale_Reptilia          0.1290  0.1533   0.8415  681  0.4004 
sex_species_classMale_Actinopterygii      0.0948  0.1901   0.4989  681  0.6180 
sex_species_classMale_Aves               -0.0352  0.0416  -0.8466  681  0.3975 
sex_species_classMale_Insecta            -0.1101  0.3352  -0.3285  681  0.7426 
sex_species_classMale_Mammalia            0.0216  0.0540   0.3994  681  0.6897 
sex_species_classMale_Reptilia            0.0132  0.2204   0.0597  681  0.9524 
                                          ci.lb    ci.ub     
sex_species_classBoth_Aves              -0.1091   0.0938     
sex_species_classBoth_Insecta           -0.6160   0.3271     
sex_species_classBoth_Mammalia          -0.3033  -0.0495  ** 
sex_species_classBoth_Reptilia          -0.3727   0.7033     
sex_species_classFemale_Actinopterygii  -0.4409   0.3224     
sex_species_classFemale_Arachnida       -0.4719   0.5539     
sex_species_classFemale_Aves            -0.1136   0.0484     
sex_species_classFemale_Insecta         -0.0981   0.5946     
sex_species_classFemale_Mammalia        -0.1102   0.0841     
sex_species_classFemale_Reptilia        -0.1720   0.4299     
sex_species_classMale_Actinopterygii    -0.2784   0.4680     
sex_species_classMale_Aves              -0.1169   0.0464     
sex_species_classMale_Insecta           -0.7684   0.5481     
sex_species_classMale_Mammalia          -0.0845   0.1277     
sex_species_classMale_Reptilia          -0.4197   0.4460     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_species_class)
   R2_marginal R2_conditional 
     0.0283798      0.3704015 
Code
orchard_plot(mod_sex_species_class, mod = "sex_species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S21. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_species_class <- all_models(mod_sex_species_class,  mod = "sex_species_class", type = "normal")

# saving tab_sex_species_class as RDS
saveRDS(tab_sex_species_class, here("Rdata", "tab_sex_species_class.RDS"))

Table S16. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and taxonomic classes with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_species_class <-readRDS(here("Rdata", "tab_sex_species_class.RDS"))

tab_sex_species_class
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_Aves -0.001 -0.116 0.114 0.987 -0.777 0.775
B_Insecta -0.146 -0.689 0.397 0.597 -1.086 0.794
B_Mammalia -0.180 -0.312 -0.047 0.008 -0.959 0.599
B_Reptilia 0.158 -0.381 0.696 0.565 -0.780 1.095
F_Actinopterygii -0.069 -0.477 0.339 0.741 -0.938 0.800
F_Arachnida 0.040 -0.514 0.594 0.887 -0.906 0.987
F_Aves -0.036 -0.124 0.051 0.417 -0.809 0.736
F_Insecta 0.240 -0.128 0.607 0.201 -0.611 1.091
F_Mammalia -0.019 -0.124 0.087 0.731 -0.793 0.756
F_Reptilia 0.071 -0.238 0.380 0.651 -0.756 0.898
M_Actinopterygii 0.068 -0.339 0.475 0.742 -0.800 0.937
M_Aves -0.035 -0.123 0.054 0.441 -0.807 0.738
M_Insecta -0.143 -0.888 0.602 0.707 -1.212 0.927
M_Mammalia -0.016 -0.129 0.097 0.778 -0.792 0.760
M_Reptilia -0.024 -0.481 0.433 0.918 -0.917 0.869
B_Aves-B_Insecta -0.145 -0.700 0.409 0.607 -1.092 0.801
B_Aves-B_Mammalia -0.179 -0.346 -0.012 0.036 -0.964 0.607
B_Aves-B_Reptilia 0.159 -0.388 0.706 0.569 -0.784 1.101
B_Aves-F_Actinopterygii -0.068 -0.489 0.353 0.752 -0.943 0.807
B_Aves-F_Arachnida 0.041 -0.523 0.605 0.886 -0.912 0.994
B_Aves-F_Aves -0.035 -0.165 0.094 0.594 -0.814 0.743
B_Aves-F_Insecta 0.241 -0.141 0.623 0.217 -0.617 1.098
B_Aves-F_Mammalia -0.018 -0.165 0.130 0.815 -0.799 0.764
B_Aves-F_Reptilia 0.072 -0.253 0.396 0.663 -0.761 0.905
B_Aves-M_Actinopterygii 0.069 -0.351 0.490 0.746 -0.806 0.944
B_Aves-M_Aves -0.034 -0.163 0.096 0.611 -0.812 0.745
B_Aves-M_Insecta -0.142 -0.894 0.611 0.712 -1.217 0.933
B_Aves-M_Mammalia -0.015 -0.169 0.138 0.845 -0.798 0.767
B_Aves-M_Reptilia -0.023 -0.491 0.445 0.923 -0.922 0.876
B_Insecta-B_Mammalia -0.033 -0.591 0.524 0.906 -0.982 0.915
B_Insecta-B_Reptilia 0.304 -0.459 1.067 0.435 -0.779 1.386
B_Insecta-F_Actinopterygii 0.077 -0.601 0.756 0.823 -0.947 1.102
B_Insecta-F_Arachnida 0.186 -0.589 0.962 0.637 -0.904 1.277
B_Insecta-F_Aves 0.110 -0.439 0.659 0.694 -0.834 1.054
B_Insecta-F_Insecta 0.386 -0.269 1.041 0.248 -0.623 1.395
B_Insecta-F_Mammalia 0.128 -0.424 0.680 0.650 -0.818 1.073
B_Insecta-F_Reptilia 0.217 -0.406 0.841 0.494 -0.771 1.206
B_Insecta-M_Actinopterygii 0.215 -0.463 0.892 0.534 -0.809 1.239
B_Insecta-M_Aves 0.112 -0.438 0.661 0.690 -0.832 1.056
B_Insecta-M_Insecta 0.004 -0.918 0.925 0.994 -1.196 1.203
B_Insecta-M_Mammalia 0.130 -0.423 0.684 0.645 -0.816 1.076
B_Insecta-M_Reptilia 0.122 -0.587 0.831 0.735 -0.922 1.167
B_Mammalia-B_Reptilia 0.337 -0.210 0.885 0.227 -0.605 1.280
B_Mammalia-F_Actinopterygii 0.111 -0.312 0.534 0.607 -0.765 0.987
B_Mammalia-F_Arachnida 0.220 -0.347 0.786 0.446 -0.734 1.174
B_Mammalia-F_Aves 0.144 -0.004 0.291 0.056 -0.638 0.925
B_Mammalia-F_Insecta 0.419 0.035 0.804 0.033 -0.439 1.278
B_Mammalia-F_Mammalia 0.161 0.020 0.303 0.026 -0.619 0.942
B_Mammalia-F_Reptilia 0.251 -0.075 0.577 0.132 -0.583 1.085
B_Mammalia-M_Actinopterygii 0.248 -0.175 0.671 0.250 -0.628 1.124
B_Mammalia-M_Aves 0.145 -0.002 0.293 0.054 -0.636 0.927
B_Mammalia-M_Insecta 0.037 -0.717 0.791 0.923 -1.039 1.113
B_Mammalia-M_Mammalia 0.163 0.015 0.312 0.031 -0.618 0.945
B_Mammalia-M_Reptilia 0.156 -0.314 0.625 0.515 -0.744 1.055
B_Reptilia-F_Actinopterygii -0.226 -0.897 0.444 0.507 -1.245 0.793
B_Reptilia-F_Arachnida -0.117 -0.887 0.652 0.764 -1.204 0.969
B_Reptilia-F_Aves -0.194 -0.735 0.347 0.482 -1.133 0.745
B_Reptilia-F_Insecta 0.082 -0.565 0.729 0.803 -0.922 1.086
B_Reptilia-F_Mammalia -0.176 -0.718 0.366 0.524 -1.116 0.764
B_Reptilia-F_Reptilia -0.087 -0.653 0.479 0.764 -1.040 0.867
B_Reptilia-M_Actinopterygii -0.089 -0.760 0.581 0.794 -1.109 0.930
B_Reptilia-M_Aves -0.192 -0.733 0.349 0.486 -1.131 0.747
B_Reptilia-M_Insecta -0.300 -1.217 0.616 0.520 -1.496 0.895
B_Reptilia-M_Mammalia -0.174 -0.718 0.370 0.530 -1.115 0.767
B_Reptilia-M_Reptilia -0.182 -0.860 0.497 0.599 -1.206 0.843
F_Actinopterygii-F_Arachnida 0.109 -0.577 0.795 0.755 -0.920 1.138
F_Actinopterygii-F_Aves 0.033 -0.381 0.446 0.877 -0.839 0.904
F_Actinopterygii-F_Insecta 0.308 -0.237 0.854 0.267 -0.633 1.250
F_Actinopterygii-F_Mammalia 0.050 -0.366 0.466 0.812 -0.823 0.923
F_Actinopterygii-F_Reptilia 0.140 -0.366 0.646 0.588 -0.779 1.059
F_Actinopterygii-M_Actinopterygii 0.137 -0.324 0.598 0.559 -0.758 1.032
F_Actinopterygii-M_Aves 0.034 -0.379 0.448 0.871 -0.838 0.906
F_Actinopterygii-M_Insecta -0.074 -0.921 0.773 0.864 -1.217 1.069
F_Actinopterygii-M_Mammalia 0.053 -0.366 0.471 0.805 -0.821 0.927
F_Actinopterygii-M_Reptilia 0.045 -0.564 0.653 0.885 -0.935 1.024
F_Arachnida-F_Aves -0.076 -0.635 0.483 0.789 -1.026 0.873
F_Arachnida-F_Insecta 0.200 -0.463 0.862 0.555 -0.814 1.213
F_Arachnida-F_Mammalia -0.059 -0.620 0.502 0.837 -1.009 0.892
F_Arachnida-F_Reptilia 0.031 -0.600 0.662 0.924 -0.963 1.024
F_Arachnida-M_Actinopterygii 0.028 -0.657 0.714 0.936 -1.001 1.057
F_Arachnida-M_Aves -0.075 -0.634 0.484 0.793 -1.024 0.875
F_Arachnida-M_Insecta -0.183 -1.110 0.744 0.699 -1.386 1.021
F_Arachnida-M_Mammalia -0.056 -0.619 0.506 0.844 -1.008 0.895
F_Arachnida-M_Reptilia -0.064 -0.780 0.651 0.860 -1.114 0.985
F_Aves-F_Insecta 0.276 -0.098 0.650 0.148 -0.578 1.130
F_Aves-F_Mammalia 0.018 -0.108 0.143 0.782 -0.760 0.795
F_Aves-F_Reptilia 0.107 -0.208 0.422 0.504 -0.722 0.937
F_Aves-M_Actinopterygii 0.104 -0.309 0.518 0.620 -0.767 0.976
F_Aves-M_Aves 0.002 -0.074 0.077 0.967 -0.770 0.773
F_Aves-M_Insecta -0.107 -0.855 0.642 0.780 -1.179 0.965
F_Aves-M_Mammalia 0.020 -0.112 0.152 0.768 -0.759 0.799
F_Aves-M_Reptilia 0.012 -0.449 0.474 0.959 -0.884 0.908
F_Insecta-F_Mammalia -0.258 -0.635 0.118 0.179 -1.113 0.597
F_Insecta-F_Reptilia -0.169 -0.643 0.305 0.485 -1.071 0.733
F_Insecta-M_Actinopterygii -0.171 -0.717 0.374 0.537 -1.113 0.770
F_Insecta-M_Aves -0.274 -0.648 0.100 0.151 -1.128 0.580
F_Insecta-M_Insecta -0.382 -1.165 0.400 0.338 -1.479 0.714
F_Insecta-M_Mammalia -0.256 -0.635 0.123 0.185 -1.112 0.600
F_Insecta-M_Reptilia -0.264 -0.846 0.318 0.374 -1.227 0.700
F_Mammalia-F_Reptilia 0.090 -0.228 0.407 0.580 -0.741 0.920
F_Mammalia-M_Actinopterygii 0.087 -0.329 0.503 0.682 -0.786 0.960
F_Mammalia-M_Aves -0.016 -0.141 0.109 0.802 -0.794 0.762
F_Mammalia-M_Insecta -0.124 -0.874 0.626 0.745 -1.197 0.949
F_Mammalia-M_Mammalia 0.002 -0.099 0.103 0.965 -0.772 0.776
F_Mammalia-M_Reptilia -0.006 -0.469 0.458 0.981 -0.902 0.891
F_Reptilia-M_Actinopterygii -0.003 -0.509 0.503 0.992 -0.922 0.917
F_Reptilia-M_Aves -0.106 -0.420 0.209 0.511 -0.935 0.724
F_Reptilia-M_Insecta -0.214 -1.017 0.590 0.602 -1.325 0.898
F_Reptilia-M_Mammalia -0.087 -0.407 0.233 0.593 -0.919 0.744
F_Reptilia-M_Reptilia -0.095 -0.542 0.352 0.676 -0.983 0.793
M_Actinopterygii-M_Aves -0.103 -0.516 0.311 0.625 -0.975 0.769
M_Actinopterygii-M_Insecta -0.211 -1.058 0.636 0.625 -1.354 0.932
M_Actinopterygii-M_Mammalia -0.085 -0.503 0.334 0.691 -0.959 0.789
M_Actinopterygii-M_Reptilia -0.092 -0.701 0.516 0.766 -1.072 0.887
M_Aves-M_Insecta -0.108 -0.857 0.640 0.777 -1.180 0.964
M_Aves-M_Mammalia 0.018 -0.114 0.151 0.786 -0.761 0.797
M_Aves-M_Reptilia 0.010 -0.451 0.472 0.965 -0.885 0.906
M_Insecta-M_Mammalia 0.126 -0.625 0.878 0.741 -0.948 1.201
M_Insecta-M_Reptilia 0.119 -0.753 0.990 0.789 -1.043 1.280
M_Mammalia-M_Reptilia -0.008 -0.473 0.457 0.974 -0.905 0.890
Code
dat$sex_whose_fitness <- as.factor(paste(dat$sex, dat$whose_fitness ,sep = "_"))

mod_sex_whose_fitness <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_whose_fitness - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_whose_fitness)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-241.1304   482.2607   500.2607   541.0909   500.5254   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0740  0.2721    696     no         effectID 
sigma^2.2  0.0205  0.1433    206     no          paperID 
sigma^2.3  0.0199  0.1410    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 690) = 15435.4374, p-val < .0001

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 690) = 1.1201, p-val = 0.3487

Model Results:

                                    estimate      se     tval   df    pval 
sex_whose_fitnessBoth_descendant      0.0341  0.1090   0.3130  690  0.7544 
sex_whose_fitnessBoth_focal          -0.0809  0.0433  -1.8666  690  0.0624 
sex_whose_fitnessFemale_descendant   -0.0108  0.0503  -0.2139  690  0.8307 
sex_whose_fitnessFemale_focal        -0.0067  0.0353  -0.1891  690  0.8501 
sex_whose_fitnessMale_descendant      0.0639  0.0568   1.1247  690  0.2611 
sex_whose_fitnessMale_focal          -0.0181  0.0364  -0.4958  690  0.6202 
                                      ci.lb   ci.ub    
sex_whose_fitnessBoth_descendant    -0.1799  0.2481    
sex_whose_fitnessBoth_focal         -0.1660  0.0042  . 
sex_whose_fitnessFemale_descendant  -0.1096  0.0881    
sex_whose_fitnessFemale_focal       -0.0759  0.0625    
sex_whose_fitnessMale_descendant    -0.0477  0.1755    
sex_whose_fitnessMale_focal         -0.0896  0.0535    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_whose_fitness)
   R2_marginal R2_conditional 
   0.009821673    0.359580965 
Code
orchard_plot(mod_sex_whose_fitness, mod = "sex_whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S22. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants. We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_whose_fitness <- all_models(mod_sex_whose_fitness,  mod = "sex_whose_fitness", type = "normal")

# saving tab_sex_whose_fitness as RDS
saveRDS(tab_sex_whose_fitness, here("Rdata", "tab_sex_whose_fitness.RDS"))

Table S17. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and focal individuals (non-dispersers and dispersers) and their descendants with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_whose_fitness <-readRDS(here("Rdata", "tab_sex_whose_fitness.RDS"))

tab_sex_whose_fitness
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_individual -0.078 -0.168 0.011 0.086 -0.849 0.693
B_offspring 0.020 -0.276 0.316 0.894 -0.801 0.841
F_individual -0.011 -0.084 0.062 0.760 -0.781 0.758
F_offspring -0.021 -0.129 0.087 0.704 -0.794 0.753
M_individual -0.034 -0.109 0.041 0.376 -0.803 0.736
M_offspring 0.061 -0.063 0.185 0.333 -0.715 0.837
B_individual-B_offspring 0.098 -0.198 0.395 0.515 -0.723 0.920
B_individual-F_individual 0.067 -0.028 0.162 0.168 -0.705 0.839
B_individual-F_offspring 0.057 -0.067 0.182 0.365 -0.719 0.833
B_individual-M_individual 0.045 -0.053 0.142 0.370 -0.728 0.817
B_individual-M_offspring 0.140 0.002 0.277 0.047 -0.639 0.918
B_offspring-F_individual -0.031 -0.329 0.266 0.836 -0.853 0.790
B_offspring-F_offspring -0.041 -0.349 0.268 0.795 -0.867 0.785
B_offspring-M_individual -0.054 -0.352 0.244 0.723 -0.876 0.768
B_offspring-M_offspring 0.041 -0.273 0.355 0.797 -0.787 0.869
F_individual-F_offspring -0.010 -0.112 0.093 0.855 -0.782 0.763
F_individual-M_individual -0.022 -0.087 0.043 0.499 -0.791 0.746
F_individual-M_offspring 0.073 -0.046 0.192 0.232 -0.703 0.848
F_offspring-M_individual -0.013 -0.118 0.092 0.810 -0.786 0.760
F_offspring-M_offspring 0.082 -0.051 0.215 0.227 -0.695 0.860
M_individual-M_offspring 0.095 -0.025 0.215 0.121 -0.680 0.870
Code
dat$sex_fitness_higher_level <- as.factor(paste(dat$sex, dat$fitness_higher_level ,sep = "_"))

mod_sex_fitness <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_fitness_higher_level - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_fitness)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-240.7301   481.4602   503.4602   553.3319   503.8507   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0738  0.2717    696     no         effectID 
sigma^2.2  0.0246  0.1569    206     no          paperID 
sigma^2.3  0.0149  0.1220    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 688) = 15355.3220, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 688) = 1.0605, p-val = 0.3891

Model Results:

                                                 estimate      se     tval   df 
sex_fitness_higher_levelBoth_reproduction         -0.0212  0.0631  -0.3357  688 
sex_fitness_higher_levelBoth_survival             -0.0917  0.0480  -1.9089  688 
sex_fitness_higher_levelFemale_indirect fitness   -0.2306  0.2135  -1.0801  688 
sex_fitness_higher_levelFemale_reproduction        0.0071  0.0368   0.1927  688 
sex_fitness_higher_levelFemale_survival           -0.0291  0.0402  -0.7230  688 
sex_fitness_higher_levelMale_indirect fitness     -0.2278  0.1748  -1.3031  688 
sex_fitness_higher_levelMale_reproduction         -0.0171  0.0399  -0.4287  688 
sex_fitness_higher_levelMale_survival              0.0163  0.0418   0.3886  688 
                                                   pval    ci.lb   ci.ub    
sex_fitness_higher_levelBoth_reproduction        0.7372  -0.1452  0.1028    
sex_fitness_higher_levelBoth_survival            0.0567  -0.1860  0.0026  . 
sex_fitness_higher_levelFemale_indirect fitness  0.2805  -0.6498  0.1886    
sex_fitness_higher_levelFemale_reproduction      0.8473  -0.0651  0.0793    
sex_fitness_higher_levelFemale_survival          0.4700  -0.1080  0.0499    
sex_fitness_higher_levelMale_indirect fitness    0.1930  -0.5710  0.1154    
sex_fitness_higher_levelMale_reproduction        0.6682  -0.0954  0.0612    
sex_fitness_higher_levelMale_survival            0.6977  -0.0659  0.0984    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_fitness)
   R2_marginal R2_conditional 
    0.01181556     0.35629363 
Code
orchard_plot(mod_sex_fitness, mod = "sex_fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Figure S23. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_fitness <- all_models(mod_sex_fitness,  mod = "sex_fitness_higher_level", type = "normal")

# saving tab_sex_fitness as RDS
saveRDS(tab_sex_fitness, here("Rdata", "tab_sex_fitness.RDS"))

Table S18. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and different fitness components (survival, reproduction, and indirect fitness) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_fitness <-readRDS(here("Rdata", "tab_sex_fitness.RDS"))

tab_sex_fitness
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_reproduction -0.023 -0.164 0.118 0.749 -0.800 0.754
B_survival -0.095 -0.196 0.006 0.066 -0.865 0.676
F_indirect fitness -0.231 -0.732 0.270 0.365 -1.144 0.682
F_reproduction -0.004 -0.081 0.073 0.922 -0.771 0.764
F_survival -0.026 -0.112 0.061 0.558 -0.794 0.743
M_indirect fitness -0.230 -0.638 0.178 0.269 -1.096 0.636
M_reproduction -0.028 -0.112 0.056 0.513 -0.796 0.740
M_survival -0.001 -0.089 0.088 0.988 -0.769 0.768
B_reproduction-B_survival -0.072 -0.231 0.088 0.378 -0.852 0.708
B_reproduction-F_indirect fitness -0.208 -0.726 0.310 0.431 -1.131 0.715
B_reproduction-F_reproduction 0.019 -0.132 0.170 0.803 -0.759 0.798
B_reproduction-F_survival -0.003 -0.159 0.153 0.972 -0.782 0.777
B_reproduction-M_indirect fitness -0.207 -0.636 0.223 0.345 -1.083 0.669
B_reproduction-M_reproduction -0.005 -0.159 0.149 0.949 -0.784 0.774
B_reproduction-M_survival 0.022 -0.135 0.180 0.781 -0.757 0.802
B_survival-F_indirect fitness -0.136 -0.644 0.371 0.598 -1.053 0.781
B_survival-F_reproduction 0.091 -0.017 0.199 0.098 -0.680 0.862
B_survival-F_survival 0.069 -0.045 0.183 0.233 -0.703 0.841
B_survival-M_indirect fitness -0.135 -0.552 0.282 0.525 -1.005 0.735
B_survival-M_reproduction 0.067 -0.046 0.180 0.247 -0.705 0.839
B_survival-M_survival 0.094 -0.023 0.211 0.114 -0.678 0.867
F_indirect fitness-F_reproduction 0.227 -0.275 0.729 0.375 -0.687 1.141
F_indirect fitness-F_survival 0.205 -0.297 0.708 0.423 -0.709 1.119
F_indirect fitness-M_indirect fitness 0.001 -0.602 0.604 0.997 -0.972 0.974
F_indirect fitness-M_reproduction 0.203 -0.300 0.706 0.428 -0.711 1.118
F_indirect fitness-M_survival 0.230 -0.272 0.733 0.369 -0.684 1.145
F_reproduction-F_survival -0.022 -0.104 0.060 0.600 -0.790 0.746
F_reproduction-M_indirect fitness -0.226 -0.636 0.184 0.279 -1.093 0.641
F_reproduction-M_reproduction -0.024 -0.102 0.054 0.542 -0.792 0.743
F_reproduction-M_survival 0.003 -0.085 0.091 0.944 -0.766 0.772
F_survival-M_indirect fitness -0.204 -0.615 0.207 0.330 -1.071 0.663
F_survival-M_reproduction -0.002 -0.093 0.089 0.962 -0.771 0.767
F_survival-M_survival 0.025 -0.064 0.114 0.579 -0.744 0.794
M_indirect fitness-M_reproduction 0.202 -0.208 0.612 0.334 -0.665 1.069
M_indirect fitness-M_survival 0.229 -0.182 0.640 0.274 -0.638 1.096
M_reproduction-M_survival 0.027 -0.067 0.122 0.570 -0.742 0.797
Code
dat$sex_dispersal_type <- as.factor(paste(dat$sex, dat$dispersal_type ,sep = "_"))

mod_sex_dispersal_type <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_dispersal_type - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex_dispersal_type)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-241.4964   482.9928   506.9928   561.3808   507.4557   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0745  0.2729    696     no         effectID 
sigma^2.2  0.0304  0.1744    206     no          paperID 
sigma^2.3  0.0087  0.0935    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 687) = 15251.0242, p-val < .0001

Test of Moderators (coefficients 1:9):
F(df1 = 9, df2 = 687) = 0.5900, p-val = 0.8059

Model Results:

                                   estimate      se     tval   df    pval 
sex_dispersal_typeBoth_Both         -0.1189  0.1211  -0.9817  687  0.3266 
sex_dispersal_typeBoth_breeding      0.0131  0.0846   0.1544  687  0.8773 
sex_dispersal_typeBoth_natal        -0.0838  0.0475  -1.7628  687  0.0784 
sex_dispersal_typeFemale_Both       -0.0864  0.1219  -0.7088  687  0.4787 
sex_dispersal_typeFemale_breeding   -0.0071  0.0498  -0.1431  687  0.8863 
sex_dispersal_typeFemale_natal      -0.0019  0.0366  -0.0518  687  0.9587 
sex_dispersal_typeMale_Both         -0.0056  0.1386  -0.0401  687  0.9680 
sex_dispersal_typeMale_breeding     -0.0257  0.0533  -0.4817  687  0.6302 
sex_dispersal_typeMale_natal         0.0083  0.0388   0.2145  687  0.8302 
                                     ci.lb   ci.ub    
sex_dispersal_typeBoth_Both        -0.3566  0.1189    
sex_dispersal_typeBoth_breeding    -0.1530  0.1791    
sex_dispersal_typeBoth_natal       -0.1772  0.0095  . 
sex_dispersal_typeFemale_Both      -0.3257  0.1529    
sex_dispersal_typeFemale_breeding  -0.1049  0.0907    
sex_dispersal_typeFemale_natal     -0.0738  0.0700    
sex_dispersal_typeMale_Both        -0.2777  0.2666    
sex_dispersal_typeMale_breeding    -0.1303  0.0790    
sex_dispersal_typeMale_natal       -0.0679  0.0845    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex_dispersal_type)
   R2_marginal R2_conditional 
   0.009425299    0.350676406 
Code
orchard_plot(mod_sex_dispersal_type, mod = "sex_dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, trunk.size = 0.5, angle = 45)

Figure S24. Mean fitness effect for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types combined). We present mean effect sizes (Zr) along with 95% confidence (thicker bar) and prediction (thinner bars) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
tab_sex_dispersal_type <- all_models(mod_sex_dispersal_type,  mod = "sex_dispersal_type", type = "normal")

# saving tab_sex_dispersal_type as RDS
saveRDS(tab_sex_dispersal_type, here("Rdata", "tab_sex_dispersal_type.RDS"))

Table S19. Mean estimates for the interaction between sex (females, males, or both sexes grouped) and natal and breeding dispersal (or both dispersal types grouped) with their 95% confidence intervals (lower and upper CI) and 95% prediction intervals (lower and upper PR) are presented along with p-value (pval).

Code
tab_sex_dispersal_type <-readRDS(here("Rdata", "tab_sex_dispersal_type.RDS"))

tab_sex_dispersal_type
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B_B -0.096 -0.321 0.129 0.403 -0.892 0.699
B_breeding -0.004 -0.200 0.192 0.967 -0.792 0.784
B_natal -0.083 -0.184 0.017 0.103 -0.853 0.686
F_B -0.081 -0.301 0.138 0.468 -0.875 0.713
F_breeding -0.004 -0.110 0.102 0.940 -0.774 0.766
F_natal -0.008 -0.086 0.070 0.840 -0.775 0.759
M_B -0.108 -0.350 0.133 0.379 -0.909 0.692
M_breeding -0.022 -0.137 0.094 0.712 -0.793 0.750
M_natal -0.004 -0.086 0.077 0.922 -0.771 0.763
B_B-B_breeding 0.092 -0.203 0.387 0.540 -0.726 0.910
B_B-B_natal 0.013 -0.223 0.249 0.916 -0.786 0.811
B_B-F_B 0.015 -0.288 0.317 0.923 -0.806 0.836
B_B-F_breeding 0.092 -0.150 0.334 0.456 -0.708 0.892
B_B-F_natal 0.088 -0.142 0.318 0.453 -0.709 0.885
B_B-M_B -0.012 -0.330 0.305 0.939 -0.839 0.814
B_B-M_breeding 0.074 -0.172 0.321 0.554 -0.727 0.876
B_B-M_natal 0.092 -0.139 0.323 0.435 -0.705 0.889
B_breeding-B_natal -0.079 -0.292 0.133 0.464 -0.871 0.713
B_breeding-F_B -0.077 -0.368 0.214 0.603 -0.893 0.739
B_breeding-F_breeding 0.000 -0.214 0.214 1.000 -0.792 0.792
B_breeding-F_natal -0.004 -0.208 0.200 0.970 -0.794 0.786
B_breeding-M_B -0.104 -0.413 0.204 0.506 -0.927 0.718
B_breeding-M_breeding -0.018 -0.238 0.202 0.875 -0.812 0.776
B_breeding-M_natal 0.000 -0.206 0.206 1.000 -0.790 0.790
B_natal-F_B 0.002 -0.232 0.237 0.985 -0.796 0.800
B_natal-F_breeding 0.079 -0.051 0.209 0.232 -0.695 0.853
B_natal-F_natal 0.075 -0.032 0.183 0.169 -0.695 0.846
B_natal-M_B -0.025 -0.281 0.231 0.848 -0.830 0.780
B_natal-M_breeding 0.062 -0.078 0.202 0.386 -0.714 0.837
B_natal-M_natal 0.079 -0.031 0.190 0.158 -0.692 0.850
F_B-F_breeding 0.077 -0.160 0.314 0.523 -0.722 0.876
F_B-F_natal 0.073 -0.153 0.300 0.526 -0.723 0.869
F_B-M_B -0.027 -0.288 0.234 0.837 -0.834 0.779
F_B-M_breeding 0.059 -0.183 0.302 0.630 -0.741 0.860
F_B-M_natal 0.077 -0.150 0.304 0.506 -0.719 0.873
F_breeding-F_natal -0.004 -0.113 0.105 0.944 -0.775 0.767
F_breeding-M_B -0.104 -0.362 0.154 0.427 -0.910 0.701
F_breeding-M_breeding -0.018 -0.130 0.095 0.758 -0.789 0.754
F_breeding-M_natal 0.000 -0.112 0.112 1.000 -0.771 0.771
F_natal-M_B -0.100 -0.349 0.148 0.427 -0.903 0.702
F_natal-M_breeding -0.014 -0.134 0.107 0.823 -0.786 0.759
F_natal-M_natal 0.004 -0.068 0.076 0.916 -0.762 0.770
M_B-M_breeding 0.087 -0.176 0.350 0.517 -0.720 0.894
M_B-M_natal 0.104 -0.145 0.354 0.412 -0.698 0.907
M_breeding-M_natal 0.018 -0.106 0.142 0.780 -0.755 0.791

5 Species-level models

There are two steps, which can be done using phylo_model().

5.1 Estimate species-level effect sizes.

Code
# function to estimate species-level effect sizes
phylo_model <- function(data = dat, tree, correlation = 0.5) {
  
  # accounting for within-species correlation
  Vm <- vcalc(vi = data$vi, cluster = data$shared_group, obs = data$effectID, rho = correlation)
  mod <- rma.mv(yi = yi, 
              V = Vm,
              mod = ~ 1,
              data = data,
              random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | phylogeny),
              #R= list(phylogeny = cor_tree),
              test = "t",
              sparse = TRUE)
  
  
  # get random effects
  # dat3 <- dat2[c("phylogeny", "yi", "vi")]
  REs <- ranef(mod)$phylogeny
  
  # calculate means and CIs
  dat <- data.frame(phylogeny = rownames(REs),
                    estimate = mod$beta[1] + REs$intrcpt,
                    se = sqrt((mod$se[1])^2 + (REs$se)^2)) %>%
    mutate(lb = estimate - se * qnorm(1 - 0.05/2), 
           ub = estimate + se * qnorm(1 - 0.05/2)) %>%
    arrange(phylogeny)
  # get species class data
  species_class <- dplyr::distinct(data, phylogeny, .keep_all = TRUE) %>%
    dplyr::select(species_class, phylogeny)
  
  # extract tip labels from the tree and merge with species class
  tip.label <- data.frame(phylogeny = tree$tip.label)  # extract tip label
  species_class2 <- left_join(tip.label, species_class, by = "phylogeny")
  
  return(list(data = dat, class = species_class2))  # return processed data and species class
}

# estimate species-level effect sizes
# you can specify your own correlation coefficient via `correlation`. For illustration purpose, I used 0.5.
result <- phylo_model(data = dat, tree = tree, correlation = 0.5)

# find the significant species
result$data <- result$data %>%
  mutate(
    p_value = 2 * (1 - pnorm(abs(estimate / se))),
    significance = ifelse(p_value < 0.05, "Significant", "NonSignificant")
  )
result$data <- result$data %>%
  mutate(significance_label = ifelse(p_value < 0.05, "*", NA))

5.2 Plot tree and add species-level effect size estimates

Code
# plot tree
p1 <- ggtree(tree, layout = "rectangular", cex = 0.3)

p2 <- p1 %<+% result$class + geom_tiplab(aes(color = species_class), size = 3, align = T, offset = 0.05) + geom_tippoint(aes(color = species_class)) + guides(color = "none")

p3 <- p2 + 
  geom_facet(
    panel = "Effect size", 
    data = result$data, 
    geom = ggstance::geom_pointrangeh,
    mapping = aes(x = estimate, xmin = lb, xmax = ub, color = species_class)
  ) +
  geom_facet(
    panel = "Effect size",
    data = result$data,
    geom = geom_text,
    mapping = aes(x = lb, label = significance_label), 
    hjust = 1.5, size = 3, color = "red"              # adjust position of star *
  ) +
  geom_facet(
    panel = "Effect size",
    data = result$data, 
    geom = geom_vline,
    mapping = aes(xintercept = 0),                    
    linetype = "dashed", color = "gray"
  ) +
  theme_tree2()
Code
# adjust widths of the plot facets as necessary, if it improves visualization 
#| fig-cap: "**Figure S25.** Phylogenetic tree (left panel) and species-level effect size estimates (right panel). Significant relationships are indicated with a red asterisk. Dispersal was positive for ruffed grouse (Bonasa umbellus; Zr = 0.42, 95% CIs: 0.16 to 0.68) and negative for brown jays (Cyanocorax morio; Zr = -0.30, 95% CIs: -0.58 to -0.01), Tasmanian nativehens (Gallinula mortierii; Zr = -0.34, 95% CIs: -0.60 to X), and cougars (Puma concolor; Zr = -0.61, 95% CIs: -0.90 to -0.32)."

facet_widths(p3, c(Tree = 0.6, `Effect size` = 0.4))

:::

6 Multi-moderator-regression model

Code
# used "maximum likelihood" for model selection

#######################
# Multi-variable models
#######################

mod_full <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1 + sex +
                age_class_clean + whose_fitness + 
                fitness_higher_level +
                dispersal_type + dispersal_phase,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              method = "ML",
              test = "t",
              sparse = TRUE)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 696; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-238.0923  2849.3379   504.1845   567.8194   504.8013   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0728  0.2699    696     no         effectID 
sigma^2.2  0.0173  0.1316    206     no          paperID 
sigma^2.3  0.0225  0.1501    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 685) = 15236.1344, p-val < .0001

Test of Moderators (coefficients 2:11):
F(df1 = 10, df2 = 685) = 1.1945, p-val = 0.2910

Model Results:

                                  estimate      se     tval   df    pval 
intrcpt                            -0.2420  0.1744  -1.3872  685  0.1658 
sexFemale                           0.0461  0.0445   1.0352  685  0.3009 
sexMale                             0.0531  0.0452   1.1761  685  0.2400 
age_class_cleanjuvenile             0.0790  0.0558   1.4171  685  0.1569 
age_class_cleanmix                 -0.0691  0.0622  -1.1124  685  0.2664 
whose_fitnessfocal                 -0.0565  0.0415  -1.3603  685  0.1742 
fitness_higher_levelreproduction    0.2257  0.1445   1.5621  685  0.1187 
fitness_higher_levelsurvival        0.1878  0.1465   1.2814  685  0.2005 
dispersal_typebreeding             -0.0120  0.0894  -0.1348  685  0.8928 
dispersal_typenatal                -0.0084  0.0817  -0.1022  685  0.9186 
dispersal_phasesettlement           0.0374  0.0501   0.7461  685  0.4559 
                                    ci.lb   ci.ub    
intrcpt                           -0.5844  0.1005    
sexFemale                         -0.0413  0.1335    
sexMale                           -0.0356  0.1418    
age_class_cleanjuvenile           -0.0305  0.1886    
age_class_cleanmix                -0.1912  0.0529    
whose_fitnessfocal                -0.1381  0.0251    
fitness_higher_levelreproduction  -0.0580  0.5095    
fitness_higher_levelsurvival      -0.0999  0.4755    
dispersal_typebreeding            -0.1875  0.1634    
dispersal_typenatal               -0.1688  0.1521    
dispersal_phasesettlement         -0.0610  0.1358    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_full)*100, 2)
   R2_marginal R2_conditional 
          2.17          36.75 
Code
# multi-model selection
candidates <- dredge(mod_full, trace = 2)

# saving candidates as RDS
saveRDS(candidates, here("Rdata", "candidates.RDS"))
Code
candidates <-readRDS(here("Rdata", "candidates.RDS"))

candidates
Global model call: rma.mv(yi = yi, V = VCV, mods = ~1 + sex + age_class_clean + 
    whose_fitness + fitness_higher_level + dispersal_type + dispersal_phase, 
    random = list(~1 | effectID, ~1 | paperID, ~1 | species_cleaned), 
    data = dat, method = "ML", test = "t", sparse = TRUE)
---
Model selection table 
   (Int) age_cls_cln dsp_phs dsp_typ ftn_hgh_lvl sex whs_ftn df   logLik  AICc
2      +           +                                          6 -276.828 565.8
4      +           +       +                                  7 -275.933 566.0
34     +           +                                       +  7 -276.278 566.7
36     +           +       +                               +  8 -275.479 567.2
1      +                                                      4 -279.768 567.6
42     +           +                           +           +  9 -274.701 567.7
10     +           +                           +              8 -275.782 567.8
18     +           +                               +          8 -275.927 568.1
12     +           +       +                   +              9 -275.053 568.4
33     +                                                   +  5 -279.217 568.5
3      +                   +                                  5 -279.276 568.6
44     +           +       +                   +           + 10 -274.319 569.0
20     +           +       +                       +          9 -275.359 569.0
50     +           +                               +       +  9 -275.518 569.3
6      +           +               +                          8 -276.641 569.5
9      +                                       +              6 -278.813 569.8
35     +                   +                               +  6 -278.815 569.8
41     +                                       +           +  7 -277.874 569.9
8      +           +       +       +                          9 -275.830 569.9
17     +                                           +          6 -278.913 570.0
26     +           +                           +   +         10 -274.949 570.2
52     +           +       +                       +       + 10 -275.002 570.3
38     +           +               +                       +  9 -276.064 570.4
5      +                           +                          6 -279.306 570.7
58     +           +                           +   +       + 11 -274.187 570.8
40     +           +       +       +                       + 10 -275.349 571.0
11     +                   +                   +              7 -278.449 571.1
49     +                                           +       +  7 -278.515 571.2
28     +           +       +                   +   +         11 -274.474 571.3
46     +           +               +           +           + 11 -274.557 571.5
19     +                   +                       +          7 -278.735 571.6
14     +           +               +           +             10 -275.674 571.7
37     +                           +                       +  7 -278.774 571.7
43     +                   +                   +           +  8 -277.769 571.8
22     +           +               +               +         10 -275.770 571.9
7      +                   +       +                          7 -278.938 572.0
25     +                                       +   +          8 -278.038 572.3
60     +           +       +                   +   +       + 12 -273.929 572.3
16     +           +       +       +           +             11 -274.996 572.4
48     +           +       +       +           +           + 12 -274.220 572.9
24     +           +       +       +               +         11 -275.261 572.9
51     +                   +                       +       +  8 -278.367 572.9
54     +           +               +               +       + 11 -275.338 573.1
57     +                                       +   +       +  9 -277.439 573.1
13     +                           +           +              8 -278.487 573.2
39     +                   +       +                       +  8 -278.490 573.2
21     +                           +               +          8 -278.571 573.4
45     +                           +           +           +  9 -277.602 573.5
27     +                   +                   +   +          9 -277.898 574.1
30     +           +               +           +   +         12 -274.856 574.2
56     +           +       +       +               +       + 12 -274.881 574.2
53     +                           +               +       +  9 -278.186 574.6
15     +                   +       +           +              9 -278.197 574.7
62     +           +               +           +   +       + 13 -274.062 574.7
59     +                   +                   +   +       + 10 -277.406 575.1
23     +                   +       +               +          9 -278.444 575.2
32     +           +       +       +           +   +         13 -274.418 575.4
47     +                   +       +           +           + 10 -277.534 575.4
29     +                           +           +   +         10 -277.795 575.9
64     +           +       +       +           +   +       + 14 -273.837 576.3
55     +                   +       +               +       + 10 -278.088 576.5
61     +                           +           +   +       + 11 -277.225 576.8
31     +                   +       +           +   +         11 -277.687 577.8
63     +                   +       +           +   +       + 12 -277.207 578.9
   delta weight
2   0.00  0.122
4   0.25  0.107
34  0.94  0.076
36  1.39  0.061
1   1.81  0.049
42  1.89  0.047
10  2.00  0.045
18  2.29  0.039
12  2.59  0.033
33  2.74  0.031
3   2.86  0.029
44  3.19  0.025
20  3.21  0.024
50  3.52  0.021
6   3.72  0.019
9   3.97  0.017
35  3.97  0.017
41  4.13  0.015
8   4.15  0.015
17  4.17  0.015
26  4.45  0.013
52  4.55  0.012
38  4.62  0.012
5   4.96  0.010
58  4.99  0.010
40  5.25  0.009
11  5.28  0.009
49  5.42  0.008
28  5.56  0.008
46  5.73  0.007
19  5.86  0.007
14  5.90  0.006
37  5.93  0.006
43  5.97  0.006
22  6.09  0.006
7   6.26  0.005
25  6.51  0.005
60  6.55  0.005
16  6.61  0.004
48  7.13  0.003
24  7.14  0.003
51  7.17  0.003
54  7.29  0.003
57  7.37  0.003
13  7.41  0.003
39  7.42  0.003
21  7.58  0.003
45  7.69  0.003
27  8.28  0.002
30  8.40  0.002
56  8.45  0.002
53  8.86  0.001
15  8.88  0.001
62  8.89  0.001
59  9.36  0.001
23  9.38  0.001
32  9.60  0.001
47  9.62  0.001
29 10.14  0.001
64 10.53  0.001
55 10.72  0.001
61 11.07  0.000
31 11.99  0.000
63 13.10  0.000
Models ranked by AICc(x) 
Code
# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 2) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 2)) 

mr_averaged_aic2

Call:
model.avg(object = candidates, subset = delta < 2)

Component model call: 
rma.mv(yi = yi, V = VCV, mods = ~<7 unique rhs>, random = list(~1 | 
     effectID, ~1 | paperID, ~1 | species_cleaned), data = dat, method = ML, 
     test = t, sparse = TRUE)

Component models: 
       df  logLik   AICc delta weight
1       6 -276.83 565.78  0.00   0.24
12      7 -275.93 566.03  0.25   0.21
14      7 -276.28 566.72  0.94   0.15
124     8 -275.48 567.17  1.39   0.12
(Null)  4 -279.77 567.59  1.81   0.10
134     9 -274.70 567.67  1.89   0.09
13      8 -275.78 567.78  2.00   0.09

Term codes: 
     age_class_clean      dispersal_phase fitness_higher_level 
                   1                    2                    3 
       whose_fitness 
                   4 

Model-averaged coefficients:  
(full average) 
                                 Estimate Std. Error z value Pr(>|z|)
intrcpt                          -0.08622    0.11029   0.782    0.434
age_class_cleanjuvenile           0.06241    0.05321   1.173    0.241
age_class_cleanmix               -0.10467    0.07067   1.481    0.139
dispersal_phasesettlement         0.02130    0.04137   0.515    0.607
whose_fitnessoffspring            0.01788    0.03542   0.505    0.614
fitness_higher_levelreproduction  0.03936    0.11059   0.356    0.722
fitness_higher_levelsurvival      0.03285    0.10117   0.325    0.745
 
(conditional average) 
                                 Estimate Std. Error z value Pr(>|z|)  
intrcpt                          -0.08622    0.11029   0.782   0.4344  
age_class_cleanjuvenile           0.06911    0.05169   1.337   0.1812  
age_class_cleanmix               -0.11591    0.06502   1.783   0.0746 .
dispersal_phasesettlement         0.06429    0.04902   1.312   0.1896  
whose_fitnessoffspring            0.04929    0.04369   1.128   0.2592  
fitness_higher_levelreproduction  0.21658    0.17002   1.274   0.2027  
fitness_higher_levelsurvival      0.18075    0.17198   1.051   0.2933  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# relative importance of each predictor for all the models
importance <- sw(candidates)

importance
                     age_class_clean dispersal_phase whose_fitness
Sum of weights:      0.74            0.40            0.40         
N containing models:   32              32              32         
                     fitness_higher_level sex  dispersal_type
Sum of weights:      0.28                 0.20 0.14          
N containing models:   32                   32   32          

7 Publication bias

Code
# if each group n is available - assume n/2
dat$effectN <- ifelse(is.na(dat$n_group_1), (dat$n/2)*2/dat$n,  
                      (dat$n_group_1 * dat$n_group_2) / (dat$n_group_1 + dat$n_group_2))
  
dat$sqeffectN <- sqrt(dat$effectN)

mod_egger <- rma.mv(yi = yi, 
                    V = VCV,
                mod = ~ sqeffectN,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_egger)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-243.1309   486.2619   496.2619   518.9742   496.3491   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0747  0.2734    696     no         effectID 
sigma^2.2  0.0225  0.1500    206     no          paperID 
sigma^2.3  0.0165  0.1283    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 15584.3157, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 0.0281, p-val = 0.8668

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -0.0282  0.0407  -0.6920  694  0.4891  -0.1082  0.0518    
sqeffectN    0.0007  0.0041   0.1677  694  0.8668  -0.0074  0.0088    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_egger)*100, 2)
   R2_marginal R2_conditional 
          0.01          34.28 
Code
small <- bubble_plot(mod_egger,
                     mod = "sqeffectN",
                     group = "paperID",
                     xlab = "sqrt(Effective N)",
                     ylab = "Effect Size: Zr",
                     g = TRUE)
Code
small

Figure S26. Relationship between sample size and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# creating publication year from "reference"
dat$reference <- as.character(dat$reference)
dat$year <- with(dat, substr(reference, nchar(reference)-4, nchar(reference)))
dat$year <- as.integer(ifelse(dat$year == "2017a" | dat$year == "2017b", 2017, dat$year))
# decline effect
mod_dec <- rma.mv(yi = yi, V = vi,
                     mod = ~ year,
                     data = dat, 
                     random = list(
                       ~ 1 | effectID,
                       ~ 1 | paperID,
                       ~ 1 | species_cleaned),
                     test = "t",
                     sparse = TRUE)
summary(mod_dec)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-223.4409   446.8818   456.8818   479.5942   456.9690   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0615  0.2480    696     no         effectID 
sigma^2.2  0.0158  0.1258    206     no          paperID 
sigma^2.3  0.0222  0.1490    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 694) = 7901.4161, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 694) = 0.4782, p-val = 0.4895

Model Results:

         estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt   -2.1368  3.0590  -0.6985  694  0.4851  -8.1427  3.8692    
year       0.0011  0.0015   0.6915  694  0.4895  -0.0019  0.0040    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_dec)*100, 2)
   R2_marginal R2_conditional 
          0.17          38.32 
Code
decline <- bubble_plot(mod_dec,
                       mod = "year",
                       group = "paperID",
                       xlab = "Publication year",
                       ylab = "Effect Size: Zr",
                       g = TRUE)
Code
decline 

Figure S27. Relationship between publication year and effect size. We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# All together

mod_comb <- rma.mv(yi = yi, 
                   V = VCV,
                   mod = ~ year + sqeffectN,
                   data = dat, 
                   random = list(
                     ~ 1 | effectID,
                     ~ 1 | paperID,
                     ~ 1 | species_cleaned),
                   # ~ 1 | phylogeny),
                   # R= list(phylogeny = cor_tree),
                   test = "t",
                   sparse = TRUE)
                   
summary(mod_comb)

Multivariate Meta-Analysis Model (k = 696; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-242.6991   485.3981   497.3981   524.6443   497.5206   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0744  0.2728    696     no         effectID 
sigma^2.2  0.0199  0.1411    206     no          paperID 
sigma^2.3  0.0210  0.1448    148     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 693) = 15568.3708, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 693) = 0.1909, p-val = 0.8263

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -1.9880  3.2601  -0.6098  693  0.5422  -8.3888  4.4128    
year         0.0010  0.0016   0.6009  693  0.5481  -0.0022  0.0042    
sqeffectN    0.0003  0.0042   0.0608  693  0.9516  -0.0079  0.0084    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_comb)*100, 2)
   R2_marginal R2_conditional 
          0.13          35.53 
Code
# combined sample size and year
bubble_plot(mod_comb,
            mod = "sqeffectN",
            group = "paperID",
            xlab = "sqrt(Effective N)",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S28. Combined effect of sample size and publication year on effect size (plotting sample size only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
# combined sample size and year
bubble_plot(mod_comb,
            mod = "year",
            group = "paperID",
            xlab = "Publication year",
            ylab = "Effect Size: Zr",
            g = TRUE)

Figure S29. Combined effect of sample size and publication year on effect size (plotting publication year only). We present mean effect sizes (Zr) along with 95% confidence (thicker dotted lines) and prediction (thinner dotted lines) intervals are depicted along with individual data points (circles) scaled by precision (1/SE). We include k = number of effect sizes (number of studies).
Code
dat <- dat %>%
  mutate(leave_out = reference)

dat$leave_out <- as.factor(dat$leave_out)


LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$leave_out))) {
  temp_dat <- dat %>%
    filter(leave_out != levels(dat$leave_out)[i])
  
  VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, obs = temp_dat$effectID, rho = 0.5)
  
  LeaveOneOut_effectsize[[i]] <-  rma.mv(yi = yi,
                                         V = VCV_leaveout, 
                                         random = list(
                                           ~ 1 | effectID,
                                           ~ 1 | paperID,
                                           ~ 1 | species_cleaned),
                                         test = "t",
                                         sparse = TRUE,
                                         data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])
}

# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(model) {
  df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)
  return(df)
}

# using dplyr to form data frame
MA_oneout <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>%
  bind_rows %>%
  mutate(left_out = levels(dat$leave_out))


# telling ggplot to stop reordering factors
MA_oneout$left_out <- as.factor(MA_oneout$left_out)
MA_oneout$left_out <- factor(MA_oneout$left_out, levels = MA_oneout$left_out)

# saving the runs
saveRDS(MA_oneout, here("Rdata", "MA_oneout.RDS"))
Code
MA_oneout <- readRDS(here("Rdata", "MA_oneout.RDS"))

# plotting
leaveoneout <- ggplot(MA_oneout) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +

  geom_hline(yintercept = mod1$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
  geom_hline(yintercept = mod1$b, lty = 1, lwd = 0.75, colour = "black") + 
  geom_hline(yintercept = mod1$ci.ub,
             lty = 3, lwd = 0.75, colour = "black") + 
  geom_pointrange(aes(x = left_out, y = est,
                      ymin = lower, ymax = upper)) + 
  xlab("Study left out") + 
  ylab("Zr (effect size), 95% CI") +
  coord_flip() + 
  theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
Code
leaveoneout

Figure S30. Leave one study out (sensitivity analysis). We present mean effect sizes (Zr) along with 95% confidence intervals.
Code
dat <- dat %>%
  mutate(leave_out = species_cleaned)

dat$leave_out <- as.factor(dat$leave_out)


LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$leave_out))) {
  temp_dat <- dat %>%
    filter(leave_out != levels(dat$leave_out)[i])
  
  VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, obs = temp_dat$effectID, rho = 0.5)
  
  LeaveOneOut_effectsize[[i]] <-  rma.mv(yi = yi,
                                         V = VCV_leaveout, 
                                         random = list(
                                           ~ 1 | effectID,
                                           ~ 1 | paperID,
                                           ~ 1 | species_cleaned),
                                         test = "t",
                                         sparse = TRUE,
                                         data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])
}

# # writing function for extracting est, ci.lb, and ci.ub from all models
# est.func <- function(model) {
#   df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)
#   return(df)
# }

# using dplyr to form data frame
MA_oneout1 <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>%
  bind_rows %>%
  mutate(left_out = levels(dat$leave_out))


# telling ggplot to stop reordering factors
MA_oneout1$left_out <- as.factor(MA_oneout1$left_out)
MA_oneout1$left_out <- factor(MA_oneout1$left_out, levels = MA_oneout1$left_out)

# saving the runs
saveRDS(MA_oneout1, here("Rdata", "MA_oneout1.RDS"))
Code
MA_oneout1 <- readRDS(here("Rdata", "MA_oneout1.RDS"))

# plotting
leaveoneout1 <- ggplot(MA_oneout1) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +

  geom_hline(yintercept = mod1$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
  geom_hline(yintercept = mod1$b, lty = 1, lwd = 0.75, colour = "black") + 
  geom_hline(yintercept = mod1$ci.ub,
             lty = 3, lwd = 0.75, colour = "black") + 
  geom_pointrange(aes(x = left_out, y = est,
                      ymin = lower, ymax = upper)) + 
  xlab("Species left out") + 
  ylab("Zr (effect size), 95% CI") +
  coord_flip() + 
  theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
Code
leaveoneout1

Figure S31. Leave one species out (sensitivity analysis). We present mean effect sizes (Zr) along with 95% confidence intervals.
Code
precision_plot<-funnel(mod1, 
       yaxis="seinv",
       # = "rstudent",
       xlab = "Standarized residuals",
       ylab = "Precision (inverse of SE)",
       ylim = c(0.001, 16),
       xlim = c(-10,15),
       digits=c(0,1)
)

Code
#| fig-cap: "**Figure S32.** Relationship between standardized residuals and precision (inverse of SE)."

precision_plot
                x          y slab
1    2.257545e-01 14.1067360    1
2    2.437112e-01 11.0905365    2
3    5.043098e-01  5.3851648    3
4   -2.428510e-01  5.8309519    4
5    1.174535e-01  5.1961524    5
6    2.630991e-01  4.5825757    6
7   -3.677233e-01  9.6436508    7
8   -5.169585e-01  4.8989795    8
9   -3.204607e-02 11.2249722    9
10  -3.722928e-02  3.3166248   10
11   2.322161e-01  2.6457513   11
12   4.522393e-01  4.4721360   12
13  -2.589598e-01  4.0000000   13
14   2.397250e-02  8.3666003   14
15  -1.452570e-01  6.4031242   15
16  -1.185648e-02  7.6157731   16
17   1.149687e-01  5.7445626   17
18  -2.679751e-02  4.6904158   18
19  -4.192155e-01  4.7958315   19
20  -5.768778e-02  7.8740079   20
21  -2.436293e-01  9.3273791   21
22  -1.506821e-01  9.2195445   22
23  -7.529412e-02  9.7979590   23
24   1.799676e-02 20.2731349   24
25   2.866794e-03 18.2756669   25
26   1.118770e-02 31.7332633   26
27   1.177102e-03 29.8998328   27
28   1.933015e-02 21.1660105   28
29   4.034446e-03 18.7616630   29
30   1.508268e-03 22.2485955   30
31  -2.039953e-03 19.4164878   31
32   8.039419e-03 32.3264598   32
33  -6.308231e-03 30.7896086   33
34  -1.547571e-03 34.4093011   34
35  -5.455978e-03 32.9696830   35
36  -2.891882e-02  9.0000000   36
37  -3.467708e-02  7.7459667   37
38   1.528532e-03  4.7958315   38
39  -8.827251e-04  4.6904158   39
40   3.435118e-02  9.2736185   40
41  -3.303843e-02  8.0622577   41
42  -8.913133e-03 14.0712473   42
43  -5.973266e-03 17.3205081   43
44  -3.992050e-03 17.6351921   44
45   5.836399e-03 22.9128785   45
46   6.519605e-03 29.8496231   46
47  -8.165449e-03 16.4316767   47
48   4.005175e-03 18.2208672   48
49   1.009986e-02 17.8325545   49
50   6.899079e-03 19.4679223   50
51   1.168969e-03 23.1732605   51
52   2.358911e-03 28.5306852   52
53   4.659848e-03 24.6373700   53
54   9.263896e-03 30.4466747   54
55  -4.276371e-02 14.7309199   55
56   1.768017e-02 17.4068952   56
57   1.668827e-01 18.0554701   57
58  -1.645347e-01 17.4068952   58
59  -2.139996e-01 18.0554701   59
60   6.938791e-02 15.3948043   60
61   2.522789e-02 45.8148448   61
62   2.496223e-02 40.3236903   62
63   2.603408e-01  7.3484692   63
64   3.092008e-01  5.0990195   64
65  -5.223968e-03  7.3484692   65
66  -1.991481e-03  8.0622577   66
67   7.769751e-03  8.6602540   67
68  -3.201348e-01  5.8309519   68
69   3.162365e-01  3.3166248   69
70  -1.600706e-02  4.5825757   70
71  -2.730220e-03  6.4031242   71
72  -1.693945e-02  6.4807407   72
73  -6.458665e-02  4.2426407   73
74   5.933512e-01  5.4772256   74
75   7.766961e-01  3.0000000   75
76   1.139562e-01  6.5574385   76
77   2.539077e-01  3.3166248   77
78   1.475377e-01  9.0553851   78
79   6.290436e-02  8.5440037   79
80  -2.094663e-02  7.3484692   80
81  -1.368791e-01  7.4161985   81
82  -8.394795e-02 10.6770783   82
83  -7.341905e-02 10.3440804   83
84   1.050567e-01 10.0000000   84
85   1.475377e-01 10.0995049   85
86   6.011110e-01  2.8284271   86
87   2.121850e-01  2.0000000   87
88   2.450645e-01  2.0000000   88
89   3.468033e-01  2.0000000   89
90  -1.047231e-02  8.2462113   90
91  -4.190936e-02  7.9372539   91
92  -8.394795e-02  6.8556546   92
93   2.094663e-02  6.8556546   93
94   1.047231e-02  2.4494897   94
95   1.796982e-01  2.0000000   95
96   6.901550e-01  2.0000000   96
97   5.870160e-01  2.0000000   97
98   2.094663e-02  8.6023253   98
99  -2.094663e-02  7.8740079   99
100  9.449314e-02  6.6332496  100
101 -3.820322e-01  6.6332496  101
102  4.676350e-01  3.7416574  102
103  1.475377e-01  3.4641016  103
104  5.593781e-01  3.4641016  104
105  3.702033e-01  3.4641016  105
106  2.707666e-01  3.1622777  106
107 -5.521915e-02 20.2484567  107
108 -8.394795e-02  7.2111026  108
109 -1.156408e-01  6.3245553  109
110 -1.156408e-01  4.5825757  110
111  4.190936e-02  4.2426407  111
112 -2.013158e-01  6.3245553  112
113 -1.050567e-01  5.4772256  113
114  8.394795e-02  4.8989795  114
115 -2.121850e-01  3.7416574  115
116 -4.658253e-02  6.1644140  116
117 -2.114431e-02  2.4494897  117
118  3.103011e-02  6.8556546  118
119 -3.446403e-02  7.2801099  119
120 -8.385307e-02  9.4868330  120
121  7.907568e-02 10.8166538  121
122 -2.602882e-02  9.4868330  122
123  2.606434e-01 10.8166538  123
124 -2.833810e-03  9.8994949  124
125 -5.308086e-03 10.2956301  125
126 -1.806303e-01  9.8994949  126
127  2.283002e-01 10.2956301  127
128  1.201397e-01  8.2462113  128
129 -8.754754e-03  8.2462113  129
130 -9.295715e-03 11.4455231  130
131 -2.949621e-02 11.4017543  131
132 -2.141748e-01 10.7703296  132
133 -3.983260e-01 10.7703296  133
134  5.915985e-03 10.7703296  134
135  4.437803e-01  7.4833148  135
136 -1.912600e-02 17.1755640  136
137  2.020382e+00  6.5574385  137
138  1.607926e+00  5.8309519  138
139  9.298626e-01  3.6055513  139
140 -2.784081e-01  5.1961524  140
141 -1.368791e-01  8.3666003  141
142  9.794204e-02  4.1231056  142
143  6.207464e-01  3.3166248  143
144  4.780581e-01  3.4641016  144
145  8.268102e-03  5.8309519  145
146  4.184211e-03  5.7445626  146
147  6.486336e-02  7.2111026  147
148 -8.205951e-03  7.2111026  148
149 -1.479444e-01  5.6568542  149
150  2.636070e-01  6.7823300  150
151  1.307399e-01  3.3166248  151
152  2.236561e-01  2.2360680  152
153  6.752949e-01 14.1774469  153
154  3.395529e-01 12.6095202  154
155  6.775988e-02 14.2478068  155
156  1.596510e-02 12.3288280  156
157  1.510655e-01  9.4868330  157
158 -3.447354e-01 10.0000000  158
159  4.329564e-02 12.1655251  159
160 -1.339904e-02 10.9544512  160
161 -2.605034e-03  9.0000000  161
162 -1.128734e+00  3.8729833  162
163 -2.550519e-01  9.1651514  163
164  1.034015e-01 11.8321596  164
165  2.152292e-01  3.4641016  165
166 -4.400714e-02 11.6189500  166
167 -6.395525e-02 11.6189500  167
168  6.615116e-02  7.8740079  168
169 -8.452585e-02  3.7416574  169
170  6.455212e-01  2.8284271  170
171 -2.598603e-02  1.7320508  171
172 -5.113620e-02  5.6568542  172
173  3.897973e-02  5.4772256  173
174  0.000000e+00  5.0000000  174
175 -9.199504e-03  3.6055513  175
176 -2.666012e-02  3.7416574  176
177  3.224649e-02  5.0000000  177
178  4.870246e-02  4.2426407  178
179  6.500626e-02  3.1622777  179
180 -5.135214e-02  4.2426407  180
181  1.585963e-01  7.0000000  181
182 -1.565522e-01 15.9059737  182
183 -3.013276e-01  3.8729833  183
184 -8.104318e-01  6.3245553  184
185 -4.573882e-01  4.3588989  185
186 -3.122930e-01  6.5574385  186
187  3.142497e-02  5.9160798  187
188 -1.047231e-02  4.5825757  188
189  3.939525e-01  4.5825757  189
190  7.077710e-02  5.9160798  190
191 -1.050567e-01  7.9372539  191
192 -3.335574e-02  6.6332496  192
193 -9.235086e-02  5.3851648  193
194 -2.571502e-02 10.1488916  194
195 -2.033839e-02  8.5440037  195
196  9.754425e-02  5.3851648  196
197  1.489083e-02  6.0827625  197
198 -3.577130e-02  8.3066239  198
199 -3.246232e-01  8.2462113  199
200 -4.456788e-01  6.3245553  200
201  1.279484e-01  4.7958315  201
202 -3.879541e-01 12.1243557  202
203 -4.028627e-01  4.6904158  203
204 -5.625189e-01  4.7958315  204
205  5.228704e-03  5.0000000  205
206  4.208413e-03 12.9614814  206
207 -3.455632e-02 12.4096736  207
208 -7.132143e-01  6.9282032  208
209 -2.111286e-02  3.6055513  209
210 -2.024994e-02  6.7823300  210
211  1.363687e-01  5.6568542  211
212  3.880550e-02  9.0000000  212
213 -1.649126e-02  8.3666003  213
214  6.266840e-02  8.3666003  214
215 -1.401986e-02 11.7898261  215
216  4.794397e-04  9.3808315  216
217 -1.396280e-02 17.2626765  217
218 -2.421086e-02 19.3907194  218
219 -1.618099e-02 15.6524758  219
220  0.000000e+00  8.0622577  220
221  2.350681e-01  8.0622577  221
222 -8.745398e-02  8.2462113  222
223  6.584093e-02  8.2462113  223
224 -5.333043e-02  5.6568542  224
225  2.440283e-01  8.1240384  225
226 -1.551663e-01  8.5440037  226
227 -1.342810e-01 12.9614814  227
228  4.948935e-02 12.3693169  228
229 -6.831082e-02 15.7480157  229
230 -6.870832e-02  6.2449980  230
231 -2.477040e-02  4.2426407  231
232  5.151882e-03  4.4721360  232
233  6.711140e-01  8.1853528  233
234 -1.176687e+00  5.5677644  234
235 -2.040857e-01  7.4161985  235
236 -1.148162e+00  4.6904158  236
237 -7.719966e-01  6.7082039  237
238  7.341905e-02  6.4807407  238
239  6.290436e-02  7.0710678  239
240 -7.341905e-02  7.9372539  240
241  5.413635e-02  8.1240384  241
242  1.047231e-02  7.2801099  242
243  2.230978e-01  7.8102497  243
244  2.340568e-01  8.2462113  244
245  2.057074e-01  8.2462113  245
246 -1.126823e-01 67.2829845  246
247 -3.420233e-02  4.0000000  247
248  3.121656e-02  4.2426407  248
249 -1.489077e-02  4.7958315  249
250 -7.677339e-03  5.1961524  250
251 -3.816861e-02  6.5574385  251
252  1.323263e-01 16.8226038  252
253  0.000000e+00 16.2480768  253
254  3.904692e-01 23.4946802  254
255 -1.210592e-01 22.7596134  255
256  1.557148e-02 18.5472370  256
257 -3.300917e-02 22.7596134  257
258 -7.899807e-02 31.7175031  258
259 -1.006306e-01 27.3861279  259
260 -2.962185e-02 31.7175031  260
261 -1.672458e-02  5.1961524  261
262 -3.185291e-02  5.4772256  262
263 -5.983457e-01  3.1622777  263
264 -1.590200e+00  2.6457513  264
265 -8.536373e-01  8.0000000  265
266 -1.288741e-01  8.0000000  266
267 -4.597819e-01  4.7958315  267
268  8.377354e-02  5.0000000  268
269 -2.020820e-01  8.2462113  269
270  1.355634e-01  9.1104336  270
271  1.107606e-01  7.3484692  271
272  3.516186e-01  6.1644140  272
273  1.267819e-01  2.6457513  273
274  4.007341e-01  5.6568542  274
275  7.622700e-03 14.8323970  275
276 -1.857110e-02  3.3166248  276
277 -1.189053e-02  2.2360680  277
278  3.617709e-02 14.8323970  278
279 -1.218842e-02 13.1909060  279
280 -8.813938e-02  3.8729833  280
281  0.000000e+00  3.1622777  281
282  3.422355e-02 13.1909060  282
283  1.981345e-01 12.6885775  283
284 -4.370226e-02 17.3493516  284
285 -1.797355e-02  7.8102497  285
286  6.962706e-03  5.9160798  286
287  1.192201e-03  6.0827625  287
288 -1.590602e-01  2.6457513  288
289  4.311082e-02  5.2915026  289
290 -1.719575e-01 17.2336879  290
291 -2.112879e-01 17.2336879  291
292 -2.450457e-01 17.2336879  292
293 -3.053223e-02 16.1554944  293
294 -4.781357e-02 16.1554944  294
295 -2.919381e-02 16.1554944  295
296 -3.459441e-02 16.1554944  296
297  1.545957e-01  4.1231056  297
298  3.445548e-01  4.1231056  298
299 -2.679346e-05  6.7823300  299
300 -8.278676e-05  4.6904158  300
301 -1.702838e+00  2.2360680  301
302 -2.146716e-01 19.9499373  302
303  5.201183e-01  6.3245553  303
304  1.475810e+00  2.0000000  304
305 -5.169245e-01  2.0000000  305
306 -1.503564e+00  2.0000000  306
307  0.000000e+00  2.0000000  307
308 -5.049834e-01  0.7071068  308
309 -1.670102e-01  1.0000000  309
310  0.000000e+00  0.7071068  310
311 -1.003041e+00  1.0000000  311
312  0.000000e+00  1.0000000  312
313  0.000000e+00  1.0000000  313
314 -1.553253e+00  2.5166115  314
315 -3.632136e-01 10.8166538  315
316 -1.169911e-01  9.5916630  316
317  5.736229e-02 11.0000000  317
318 -5.981284e-02 10.7703296  318
319 -2.774061e-03  7.6811457  319
320  2.312782e-02 20.3715488  320
321 -1.854266e-01 20.3715488  321
322 -3.128137e-01 13.3416641  322
323  3.575758e-02 20.7123152  323
324  3.516576e-02 30.6431069  324
325  1.219666e-01  5.6568542  325
326  1.141158e-01  5.0990195  326
327 -2.341895e-01  5.3851648  327
328 -5.004173e-02  6.4031242  328
329 -3.428283e-01  6.3245553  329
330  5.331659e-02  6.3245553  330
331  1.314662e-01  6.3245553  331
332 -1.927197e-01  6.9041051  332
333 -2.132695e-01  6.9041051  333
334 -6.260941e-01  6.9041051  334
335 -7.149278e-01  6.9041051  335
336 -6.047526e-02 22.8910463  336
337  0.000000e+00 14.9666295  337
338  6.628107e-01 15.9059737  338
339  6.783292e-01 15.9059737  339
340  3.929165e-01 10.8166538  340
341  3.997905e-01 10.8166538  341
342  9.940676e-01 17.2336879  342
343  6.852754e-01 17.2336879  343
344  6.450918e-01 11.0000000  344
345  6.730322e-01 11.0000000  345
346 -2.646652e+00  8.8317609  346
347 -3.516632e-01 10.3440804  347
348 -4.422286e-01 10.3440804  348
349 -2.941672e-01 13.0384048  349
350 -4.578327e-03 15.9373775  350
351 -1.826024e-01 11.4017543  351
352 -1.437723e-01  9.7467943  352
353 -1.925989e-01  9.7467943  353
354 -1.745748e-01 12.5299641  354
355  2.598932e-02 15.5241747  355
356 -2.211041e-01 10.4880885  356
357 -1.424252e-01 14.2126704  357
358  3.526692e-01  5.3851648  358
359 -1.777979e-02  9.5393920  359
360 -7.187705e-04  9.5393920  360
361 -6.117307e-03  9.5393920  361
362  1.592056e-01 10.5356538  362
363 -6.017250e-04 10.5356538  363
364 -5.121174e-03 10.5356538  364
365 -3.134858e-01 10.0995049  365
366  3.199573e-02 10.0000000  366
367  1.769198e-03 10.0000000  367
368 -1.118788e-03 10.0000000  368
369  9.796960e-02 10.7703296  369
370  1.685602e-03 10.7703296  370
371 -1.065924e-03 10.7703296  371
372  1.653262e+00 10.1980390  372
373  2.646652e+00  8.4852814  373
374 -3.496485e-01  6.6332496  374
375  1.343952e+00  4.5825757  375
376 -2.646652e+00  2.8284271  376
377  1.462460e-01 18.1383571  377
378  1.162188e-02  9.0553851  378
379  6.157414e-02  8.2462113  379
380 -6.495771e-01  6.3245553  380
381  1.336082e-01  4.1231056  381
382  1.244478e+00  2.2360680  382
383 -2.527649e-01 12.5299641  383
384 -4.856020e-02 12.5299641  384
385  1.552255e-01  9.6953597  385
386  3.055843e-01  9.6953597  386
387  1.295430e-01 12.2474487  387
388 -5.930395e-02 12.2474487  388
389  4.517174e-03 10.2956301  389
390 -1.737471e-01 10.2956301  390
391  8.298602e-02 10.1488916  391
392 -5.234141e-03  6.0000000  392
393 -1.224991e-02  8.1240384  393
394  4.037179e-01 11.0000000  394
395  1.245028e-01 13.8924440  395
396  3.172745e-01 15.7480157  396
397  3.350284e-01  6.6332496  397
398  2.629967e-01  8.4261498  398
399 -1.223592e-04  5.0000000  399
400  3.092570e-03  5.0000000  400
401 -8.237875e-01  5.0000000  401
402 -5.434200e-02  5.9160798  402
403 -1.487834e-01  9.3808315  403
404  5.264227e-02  9.3808315  404
405 -7.578089e-02  9.4868330  405
406  0.000000e+00  9.4868330  406
407  1.324651e-01  9.0000000  407
408  4.840012e-02  9.0000000  408
409 -2.390871e-01  9.3808315  409
410 -3.621135e-01  9.3808315  410
411 -1.419281e-01 12.2882057  411
412 -4.750497e-02  6.0000000  412
413 -7.083210e-02  6.0000000  413
414 -6.981969e-02 13.1148770  414
415  3.795211e-02 15.4919334  415
416 -2.041259e-01 12.2474487  416
417 -2.738542e-01  9.1104336  417
418  8.544868e-02 54.5710546  418
419  5.671018e-02 48.8466990  419
420  0.000000e+00  7.7459667  420
421 -1.685145e-02 10.4880885  421
422 -5.393746e-01  7.7459667  422
423 -1.060603e-01 10.4880885  423
424 -9.898846e-01  7.7459667  424
425 -6.607228e-01 10.4880885  425
426 -2.566962e-01  4.6368092  426
427 -1.284211e-01  4.6368092  427
428  2.866274e-01  8.9442719  428
429  4.220343e-01  8.6023253  429
430  9.219314e-02 14.3874946  430
431  1.181346e-01 15.5563492  431
432 -3.128739e-01 21.2837967  432
433 -4.040298e-02 18.4661853  433
434  5.569616e-02 20.0066316  434
435 -1.046453e-02 17.3556241  435
436 -8.849302e-02  8.8317609  436
437  7.527060e-02 10.1980390  437
438 -8.679918e-03  9.6436508  438
439 -7.424981e-03 10.9544512  439
440  4.288173e-02 14.5945195  440
441  1.823140e-01 14.5945195  441
442  9.408067e-02 19.2613603  442
443  9.644150e-02 19.2613603  443
444 -6.330797e-01 10.2469508  444
445 -4.184568e-01  9.6953597  445
446  1.085125e-01  5.7445626  446
447  7.892888e-02  6.7082039  447
448 -2.334742e-02 21.2602916  448
449 -1.906140e-01 21.4941853  449
450 -1.191073e-01 21.4941853  450
451  1.020483e-01  1.7320508  451
452 -2.228607e-01 19.5448203  452
453 -1.950404e-01 19.7484177  453
454  6.196177e-02 20.6639783  454
455  1.215160e-03 20.5669638  455
456  5.872766e-02 20.9523268  456
457  1.058292e-02 20.8326667  457
458  1.123188e-01 17.2046505  458
459  7.481341e-02 16.9410743  459
460  4.860607e-02 11.5325626  460
461  9.674692e-02 12.8062485  461
462  3.112705e-01  5.7445626  462
463 -2.543599e-01  3.7416574  463
464  6.655062e-01  4.0000000  464
465  1.531290e-02  7.6811457  465
466 -7.387549e-01  5.4772256  466
467  1.740500e-01  7.2801099  467
468 -3.519079e-01 12.7671453  468
469 -2.960508e-02  8.3066239  469
470 -2.891905e-01 10.6301458  470
471 -1.896796e-01 10.3923048  471
472  1.293031e+00 12.0000000  472
473  7.403535e-01  5.4772256  473
474  2.421097e-01  5.5856960  474
475  4.702252e-01  3.8013156  475
476 -2.501492e-02 22.6053091  476
477 -4.619155e-03 24.7991935  477
478 -2.184298e-03  8.4270097  478
479 -2.364723e-01 28.9827535  479
480 -2.189586e-01 28.9827535  480
481 -4.070719e-01 28.9827535  481
482 -3.485906e-01 28.9827535  482
483  2.966495e-01 20.7123152  483
484  3.178898e-01 20.7123152  484
485  0.000000e+00 24.0416306  485
486 -8.569546e-02 24.0416306  486
487 -6.670566e-01 24.0416306  487
488 -9.106068e-02 24.0416306  488
489  1.193549e-01 13.0384048  489
490 -4.107452e-01 13.0384048  490
491  2.201255e+00 13.0384048  491
492 -5.704623e-02 13.0384048  492
493 -1.607204e-01 13.6014705  493
494 -1.626318e-01 13.6014705  494
495  1.478635e-01 11.6189500  495
496 -2.285130e-01 11.6189500  496
497  4.318917e-02 13.1529464  497
498  8.084761e-02  7.8740079  498
499  1.336821e-01  2.2360680  499
500 -6.777814e-02 17.8325545  500
501  4.891945e-02 11.8743421  501
502  1.756371e-01  5.3851648  502
503 -6.645371e-02  5.3851648  503
504 -1.803995e-01 21.9544984  504
505 -1.630099e-01 22.5166605  505
506  8.177266e-02  8.0622577  506
507  2.417408e-02 10.0498756  507
508  3.708263e-01 10.0498756  508
509  3.028810e-01 10.0498756  509
510  8.514143e-03  4.0000000  510
511 -3.121375e-02  4.1231056  511
512  8.294336e-03 15.1986842  512
513 -1.929462e-01  5.2915026  513
514 -2.125337e-01  4.3588989  514
515  3.852158e-01 12.0415946  515
516  4.357086e-01  8.5440037  516
517 -1.125932e-02 17.4355958  517
518  2.247001e-02 17.4355958  518
519 -9.128088e-03 17.4355958  519
520  1.057836e-02 17.4355958  520
521 -2.216105e-01 10.6301458  521
522 -2.284486e-01 10.6301458  522
523 -1.694410e-01 10.6301458  523
524 -1.437463e-01 11.6189500  524
525 -1.387273e-01 11.6189500  525
526 -1.612642e-01 11.6189500  526
527  1.407260e-01 25.3574447  527
528  6.895828e-02 25.3574447  528
529  3.463771e-02 25.3574447  529
530 -5.835267e-02 25.3574447  530
531 -4.972787e-02 48.4974226  531
532 -4.041147e-02 24.9599679  532
533  2.633426e-02 24.9599679  533
534 -4.369849e-03 24.9599679  534
535  4.297921e-02 24.9599679  535
536  0.000000e+00 48.4974226  536
537 -1.274849e-01 12.7279221  537
538  1.826469e-03 16.5227116  538
539 -2.642477e-01 16.5227116  539
540  1.725424e-02  3.4641016  540
541 -2.583036e-02  4.0000000  541
542  3.047749e-01  3.4641016  542
543  1.155183e-01  3.1622777  543
544 -2.115356e-03 16.5227116  544
545 -9.958176e-03 19.8746069  545
546  3.209840e-01  7.9372539  546
547 -6.653686e-02 49.8297100  547
548 -8.674014e-02 49.8397432  548
549 -3.520629e-02 49.8297100  549
550 -3.564537e-02 49.8397432  550
551 -9.846450e-03 11.0000000  551
552 -5.670687e-02 11.4017543  552
553 -4.991745e-02 11.4017543  553
554  1.658721e-02 11.4017543  554
555 -1.133863e-01 11.0453610  555
556  2.174810e-02 12.2882057  556
557  1.261697e-01 10.6770783  557
558 -9.604962e-02 10.6770783  558
559  1.136514e-03 10.6770783  559
560 -3.209264e-02 12.7279221  560
561  3.577516e-01 10.5356538  561
562  1.226545e-03 10.8627805  562
563  5.750342e-04 10.7238053  563
564 -2.451252e-01 12.4498996  564
565 -2.550714e-02 12.6095202  565
566  3.371423e-01 12.6095202  566
567 -3.383595e-01  7.4833148  567
568 -4.198831e-01  6.7823300  568
569  0.000000e+00 12.6095202  569
570  1.544381e-01 12.4498996  570
571  1.088709e-01 12.4498996  571
572  6.896537e-02 11.4891253  572
573  4.998827e-02 13.6014705  573
574  1.074805e-02 13.6014705  574
575  0.000000e+00 12.9228480  575
576 -5.466672e-02 10.7703296  576
577 -3.034590e-02 10.7703296  577
578 -1.335185e-01 10.1488916  578
579 -1.393711e-01 12.0000000  579
580 -1.755145e-01 12.0000000  580
581 -3.004416e-01 11.3578167  581
582  3.330564e-02  7.2111026  582
583  2.120167e-01 10.5356538  583
584 -5.276873e-02 17.7763888  584
585  1.315573e-01 10.5356538  585
586 -1.352162e-01 17.7763888  586
587  1.934281e-01  7.1414284  587
588 -4.272568e-02 17.4642492  588
589  1.479438e-01  7.1414284  589
590 -2.009451e-02 17.4642492  590
591 -2.990426e-01  5.3851648  591
592  1.939476e-01 27.3861279  592
593  5.929142e-01 24.5153013  593
594  1.771912e-01  2.8284271  594
595 -6.780785e-02 28.2134720  595
596 -9.902167e-02 21.2602916  596
597 -2.087050e-01 22.1359436  597
598 -3.787194e-02  9.4339811  598
599 -3.248448e-03 11.4891253  599
600  2.239207e-03  6.4031242  600
601 -3.717602e-02 25.7681975  601
602  1.715945e-01 21.1660105  602
603  1.043436e-01 12.2065556  603
604  1.093360e-01  9.5916630  604
605 -7.834137e-02 19.7989899  605
606  2.568295e-01 19.7989899  606
607 -2.267395e-02 19.7989899  607
608  7.766600e-02 12.4498996  608
609 -1.416928e-02 12.4498996  609
610 -2.251844e-01 11.1355287  610
611 -1.985923e-01 11.1355287  611
612 -2.660220e-01 41.7731971  612
613 -1.572387e-01 41.7731971  613
614 -5.340650e-02 41.7731971  614
615 -4.651090e-01 40.7676342  615
616 -3.288808e-01 40.7676342  616
617 -1.473208e-01 40.7676342  617
618  3.078085e-01  8.5440037  618
619 -2.921410e-01  8.5440037  619
620  3.297250e-01  8.5440037  620
621  1.168035e-01 66.2117814  621
622 -9.383592e-03 58.6003413  622
623 -4.332986e-02 58.2923666  623
624 -5.998452e-02 62.8012739  624
625 -9.842192e-02 55.0363516  625
626 -3.820617e-03 54.7357287  626
627  1.145493e-01 65.6962708  627
628  8.793323e-03 54.7813837  628
629 -1.693720e-02 54.4058820  629
630 -4.359594e-02 64.1092817  630
631 -9.119118e-02 52.9905652  631
632  3.057377e-02 52.6307895  632
633 -4.583225e-01  5.5677644  633
634 -1.502894e-01 15.0665192  634
635 -1.003906e-01  5.2915026  635
636 -2.336321e-01  5.3851648  636
637 -3.557587e-02 21.4709106  637
638 -2.715112e-02 21.4709106  638
639 -3.202322e-02 22.6274170  639
640 -5.563761e-02 22.6274170  640
641 -7.525130e-04 15.0000000  641
642  1.556088e-01 14.1421356  642
643  5.954508e-03 14.1421356  643
644 -1.460744e-01 47.1168760  644
645 -3.712330e-02 26.4952826  645
646  3.378140e-02 33.1360831  646
647  6.757283e-02 34.4383507  647
648  6.290436e-02 14.1067360  648
649 -5.240181e-02 13.2287566  649
650 -1.368791e-01 12.7279221  650
651  1.047231e-02 11.5758369  651
652 -1.582255e-01 13.3041347  652
653  1.047231e-02 14.8323970  653
654 -1.475377e-01 13.3790882  654
655 -1.156408e-01 13.3790882  655
656 -3.142497e-02 12.2474487  656
657 -2.094663e-02 14.0000000  657
658 -2.977186e-03 17.6068169  658
659 -1.841452e-03 18.7616630  659
660  2.648531e-02 12.0000000  660
661 -5.958172e-02 12.5698051  661
662  6.521182e-02 47.7702836  662
663 -4.365255e-02 71.2530701  663
664  4.897141e-02 13.5277493  664
665  3.184842e-02 17.7200451  665
666  1.405271e-01 12.4498996  666
667  6.175762e-02 15.0665192  667
668  1.367289e-01 13.7840488  668
669 -7.915869e-02  9.5916630  669
670  1.684283e-02  5.3851648  670
671 -8.777618e-03  8.6023253  671
672 -5.460289e-02  3.8729833  672
673 -6.774275e-03  4.4721360  673
674  1.009432e-03  6.0000000  674
675  1.001839e-03  7.9372539  675
676  0.000000e+00  8.7177979  676
677  1.880004e-01  9.1651514  677
678 -4.660494e-02  8.7177979  678
679  1.949342e-01  9.1651514  679
680  1.296472e-02  8.7177979  680
681 -2.006654e-02  9.1651514  681
682  4.098572e-01  3.1622777  682
683 -2.797457e-01  2.6457513  683
684 -3.258291e-01  6.5574385  684
685 -9.621922e-01  2.0000000  685
686 -2.090577e-01  6.5574385  686
687 -2.061022e-01  3.3166248  687
688 -1.827741e-01  5.0000000  688
689  3.779087e-02 30.7733651  689
690  7.096320e-02 29.6479342  690
691 -6.911093e-03 36.0138862  691
692 -1.791255e-02 35.0428309  692
693  7.371217e-02 34.0881211  693
694  1.416174e-01 27.5680975  694
695  5.631242e-03 28.9309523  695
696 -2.867749e-02 23.7907545  696

8 R Session Information

Code
# pander for making it look nicer
sessionInfo() %>% pander()

R version 4.4.3 (2025-02-28)

Platform: x86_64-apple-darwin20

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: cowplot(v.1.1.3), matrixcalc(v.1.0-6), orchaRd(v.2.0), here(v.1.0.1), patchwork(v.1.3.0), ggstance(v.0.3.7), ggtree(v.3.14.0), ape(v.5.8-1), metafor(v.4.8-0), numDeriv(v.2016.8-1.1), metadat(v.1.4-0), Matrix(v.1.7-2), pander(v.0.6.5), gridExtra(v.2.3), emmeans(v.1.10.7), MuMIn(v.1.48.4), kableExtra(v.1.4.0), lubridate(v.1.9.4), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.4), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): tidyselect(v.1.2.1), viridisLite(v.0.4.2), vipor(v.0.4.7), farver(v.2.1.2), latex2exp(v.0.9.6), fastmap(v.1.2.0), lazyeval(v.0.2.2), TH.data(v.1.1-3), pacman(v.0.5.1), mathjaxr(v.1.6-0), digest(v.0.6.37), estimability(v.1.5.1), timechange(v.0.3.0), lifecycle(v.1.0.4), survival(v.3.8-3), tidytree(v.0.4.6), magrittr(v.2.0.3), compiler(v.4.4.3), rlang(v.1.1.5), tools(v.4.4.3), yaml(v.2.3.10), knitr(v.1.49), labeling(v.0.4.3), htmlwidgets(v.1.6.4), xml2(v.1.3.7), aplot(v.0.2.4), multcomp(v.1.4-28), withr(v.3.0.2), grid(v.4.4.3), stats4(v.4.4.3), xtable(v.1.8-4), colorspace(v.2.1-1), scales(v.1.3.0), MASS(v.7.3-64), cli(v.3.6.4), mvtnorm(v.1.3-3), rmarkdown(v.2.29), treeio(v.1.30.0), generics(v.0.1.3), rstudioapi(v.0.17.1), tzdb(v.0.4.0), ggbeeswarm(v.0.7.2), splines(v.4.4.3), parallel(v.4.4.3), ggplotify(v.0.1.2), vctrs(v.0.6.5), yulab.utils(v.0.1.9), sandwich(v.3.1-1), jsonlite(v.1.9.0), gridGraphics(v.0.5-1), hms(v.1.1.3), beeswarm(v.0.4.0), systemfonts(v.1.2.1), glue(v.1.8.0), codetools(v.0.2-20), stringi(v.1.8.4), gtable(v.0.3.6), munsell(v.0.5.1), pillar(v.1.10.1), htmltools(v.0.5.8.1), R6(v.2.6.1), rprojroot(v.2.0.4), evaluate(v.1.0.3), lattice(v.0.22-6), ggfun(v.0.1.8), Rcpp(v.1.0.14), svglite(v.2.1.3), coda(v.0.19-4.1), nlme(v.3.1-167), mgcv(v.1.9-1), xfun(v.0.51), fs(v.1.6.5), zoo(v.1.8-13) and pkgconfig(v.2.0.3)